[type=editor,auid=000,bioid=1] [1] [auid=000,bioid=1] [auid=000,bioid=1]

1] organization=Workgroup for Infrastructure Policy (WIP), Technische
Universität Berlin, addressline=Straße des 17. Juni 135, city=10623 Berlin,
country=Germany 2] organization=Energy, Transportation, Environment
Department, German Institute for Economic Research (DIW Berlin),
addressline=Mohrenstraße 58, city=10117 Berlin, country=Germany 3]
organization=Energy and Process Systems Engineering, Department of Mechanical
and Process Engineering, ETH Zurich, addressline=Tannenstrasse 3, city=Zurich
8092, country=Switzerland

[cor1]Corresponding author.

# Stabilized Benders decomposition for energy planning under climate
uncertainty

Leonard Göke<EMAIL_ADDRESS>Felix Schmidt  Mario Kendziorski  [  [  [

###### Abstract

This paper applies Benders decomposition to two-stage stochastic problems for
energy planning under climate uncertainty, a key problem for the design of
renewable energy systems. To improve performance, we adapt various refinements
for Benders decomposition to the problem’s characteristics—a simple continuous
master-problem, and few but large sub-problems. The primary focus is on
stabilization methods, specifically comparing established bundle methods to a
quadratic trust-region approach for continuous problems.

In a quantitative case-study stabilization accelerates Benders decomposition
by a factor of four to ten. Among the investigated stabilization methods, the
trust-region approach demonstrates the best performance, followed by the level
bundle method, which is on average 33% slower. The use of inexact cuts
achieved an average reduction of computation time by 28%. Thanks to parallel
computing our refined algorithm outperforms the deterministic equivalent if
planning covers more than four weather years.

Finally, we propose additional enhancements to the algorithm. First, heuristic
methods to narrow down the solution space of the master-problem. Second,
approximations of the sub-problems to determine inexact yet valid cuts faster.

###### keywords:

OR in energy Large scale optimization Benders decomposition Stabilization
Renewable energy Climate uncertainty

##  1 Introduction

To achieve the objectives of the Paris Climate Agreement and limit global
warming, the energy system must undergo a major transformation and replace 86%
of primary energy from fossil fuels as of 2019 with renewable resources
(Ritchie Roser, 2020). At the heart of this transformation is renewable
electricity from wind and photovoltaic (PV) for two reasons: First, its
technical potential exceeds other renewables and even demand projections
(Creutzig ., 2017). Estimates for the global potential of PV alone range from
1,585 to 50,580 EJ, at least three times 2019’s primary consumption. Second,
the use of renewable electricity is not limited to the power sector, but can
also be deployed for electric heating, mobility, and the creation of synthetic
fuels to decarbonize the heat, transport, and industry sector.

###  1.1 System planning under climate uncertainty

The key tools for analyzing the transformation towards renewable electricity
are techno-economic planning models, stylized representations of macro-energy
systems describing the flow and conversion of energy. Formulated as linear
optimization problems, these models decide from a portfolio of technologies on
their expansion and operation to satisfy an exogenous final demand at minimum
costs and subject to boundary conditions, for instance, emission limits.
Typical applications include the development of long-term scenarios for the
energy system, assessment of technologies in a system context, and analysis of
energy policy (Göke ., 2023).

The fluctuating nature of PV and wind generation is a challenge for planning
models. For conventional energy systems characterized by dispatchable
generation from thermal power plants, a small number of representative time-
periods is sufficient to achieve accurate results. Capturing fluctuations of
PV and wind however requires a much higher temporal resolution, that greatly
increases model size and easily renders the linear optimization problem
computationally intractable (Göke Kendziorski, 2021). Addressing this issue,
different techniques have been proposed to reduce computational complexity
while representing renewables accurately, for instance iteratively adjusting
the representative time-periods or limiting high resolution to selected parts
of the system (Teichgraeber ., 2021; Göke, 20212).

Overall efforts to make models computationally tractable and capture
fluctuations of renewables are focused on representing a single weather year
accurately, but research finds substantial variability across different
weather years, for instance in renewable generation, that heavily impact model
results (Pfenninger Staffell, 2016; Pfenninger, 2017). Several studies analyze
how the consideration of multiple weather years affects planning power
systems. Using 40 weather years, HC. Bloomfield . (2016) investigate effects
on the British power system and conclude robust planning should consider 10
weather years at least. As a key driver of inter-annual variability, they
identify the variation of wind generation. Collins . (2018) compare the impact
of different weather years on the European power system finding that ambitious
expansion of renewables will increase the climate-driven variation of
emissions and generation costs by a factor of 5 until 2030 compared to 2015.
Ruhnau Qvist (2022) study the effect of inter-annual variability on storage
requirements for a stylized fully renewable German power system. Compared to a
single average weather year, the amount of energy to be stored more than
doubles when considering 35 years of climate data. Beyond the power sector,
Lombardi . (2022) observe great variations of heating demand across weather
years for Italy. Ohlendorf Schill (2020) specifically analyze the frequency of
low-wind power events in Germany threatening system adequacy, again concluding
planning models should consider more weather years.

Heuristic approaches to consider climate uncertainty are limited. Solving a
single-year planning model for various weather years can indicate the level of
variability, but not provide an adequate and cost-efficient solution across
all years. Using representative time-periods but basing them on multiple
weather years cannot fully capture fluctuations of renewables (Hilbers .,
2019).

A more robust and analytical approach to climate variability is stochastic
programming. Formulated as two-stage problems, planning models decide on the
expansion of technologies in the first stage and on their operation in the
second. Operation in the second stage is subject to uncertainty represented by
scenarios that reflect weather-dependency and are based on distinct years of
consistent climate data. On the downside, multiple scenarios, representing one
weather year each, greatly increase model-size and aggravate the problem of
computational tractability described above. For instance, to enable two-stage
planning the EMPIRE model uses representative time-periods and omits seasonal
storage, a key component for renewable energy systems (Backe ., 2022).

###  1.2 Benders decomposition for two-stage stochastic problems

Benders decomposition (BD), first introduced in Benders (1962), is a
decomposition technique to solve optimization problems. Van Slyke Wets (1969)
first applied a variation of BD, termed the L-shaped method, to solve two-
stage stochastic problems and potentially improve their computational
tractability. Applied to energy planning with multiple scenarios, or weather
years, the problem is decomposed into a master-problem (MP) addressing
technology expansion, and a number of mutually independent sub-problems (SPs)
for operation. Afterward, the MP and SPs are solved repeatedly to generate
constraints, so-called Benders cuts, that are added to the MP until the
algorithm converges (Conejo ., 2006).

In several applications of BD for energy planning models, each SP corresponds
to a short time-span of the year, in some models as short as an hour (Lohmann
Rebennack, 2017; Brandenberg Stursberg, 2021; Skar ., 2014). However, such
decomposition is not applicable when modeling high renewable shares and
including seasonal storage, which creates dependencies across the entire year.
Thus, an accurate representation of these dependencies requires that each SP
covers operation for a specific weather year or scenario instead. Considering
the high temporal detail of operation, this results in a simple MP and several
large SPs. Furthermore, MP and SPs can contain both discrete and continuous
variables, but in the MP at least renewable expansion is continuous
considering the size of a single wind turbine or PV panel is small from a
macro perspective.

For many problems, the original BD converges slowly and is not competitive to
off-the-shelf solvers for the deterministic equivalent. However, an extensive
branch of research proposes different enhancements to improve the original BD.
In a review, Rahmaniani . (2017) group them into the following three most
important groups111We omitted the fourth group ”Decomposition strategies”
since none of the few approaches in this group is relevant for us.:

  * •

Solution procedure: The majority of enhancements in this area reduce the
computation time of the MP, like removing non-binding cuts or not solving to
optimality in early iterations. The at least in our case more relevant SPs can
deploy distributed parallel computing or sensibly reduce the number of
scenarios in early iterations. Furthermore, not solving SPs to optimality can
produce inexact but valid cuts (Zakeri ., 2000).

  * •

Solution generation: Again, most enhancements regarding solution generation
aim to accelerate the MP, for instance by solving a linear relaxation instead
of a mixed-integer problem. Other adjustments to the MP improve the
convergence rate, rendering them sensible even if the MP is small compared to
the SPs. These include:

    * –

Multi-cut reformulations that use separate cuts for each SP instead of a
single aggregated cut improving convergence but making the MP harder.

    * –

Problem-specific heuristics that provide an initial feasible solution for BD,
obtain an upper bound on the objective value, or exclude infeasible solutions
at an early stage. To achieve the latter, additional constraints, known as
valid inequalities, which are guaranteed to hold true in the optimum, are
added to the MP narrowing down its solution space (Cordeau ., 2006). For
instance, Wang . (2013) add the second-stage constraints of one scenario to
the MP and improve convergence as a result.

    * –

Regularization restricting the MP to solutions close to a reference point,
usually the current best solution, to avoid heavy oscillation (Ruszczyński,
2003). These methods are not limited to BD or stochastic programming and
rather have their origins in general convex optimization. For instance,
Linderoth Wright (2003) apply the $l_{\infty}$-norm to limit the maximum
difference across all variables between a new and the current best solution of
the MP; a method that has been previously introduced in a general optimization
context in Marsten . (1975). Level and proximal bundle methods originate from
nonsmooth convex optimization and add a quadratic term to the objective
function for stabilization (Ruszczyński, 1986; Lemaréchal ., 1981). A series
of publications specifically applies these methods to mitigate the oscillation
of BD (e.g. van Ackooij de Oliveira,, 2015; Zverovich ., 2012).

  * •

Cut generation: Magnanti Wong (1981) note that degenerate SPs, e.g. SPs
without unique solutions, can generate a set of different cuts, all affecting
convergence differently. To select strong cuts improving convergence, they
propose to solve a modified version of the dual SP after the original SP. As a
result, the method is most efficient, if the MP is difficult and the SPs are
small.

In conclusion, many of the established enhancements aim at problems with a
difficult MP and small SPs, typically combinatorial problems with discrete
complicating variables. Rahmaniani . (2017) correspondingly state ”it has
often been reported that more than 90% of the total execution […] time is
spent on solving the MP”. As a result, applicability to the two-stage planning
problem investigated in this paper with a simple MP, but large SPs is limited.
In contrast to typical applications of BD, it is not combinatorial complexity
but the mere size of the individual subproblems that make our case
computationally challenging.

###  1.3 Contribution

In this paper, we apply BD to solve a two-stage stochastic problem with
discrete scenarios used for planning renewable macro-energy systems under
climate uncertainty. The existing literature neither addresses this problem
specifically nor generally covers BD for continuous linear problems (LPs) that
are challenging due to the size of the individual subproblems.

The main contribution of the paper is to adapt and apply existing refinements
of BD to the given problem. Overall, our contribution is rather practical than
theoretical and has a focus on efficiently solving the problem at hand. The
examined refinements can be divided into two groups.

  1. 1.

First, common improvements for BD including inexact cuts, valid inequalities,
and parallelization. Many other state-of-the-art refinements are not examined
since they do not suit the specific structure of the problem.

  2. 2.

Second, quadratic stabilization, or regularization, methods that originate
from nonsmooth convex optimization and generally address the instability of
cutting plane methods. A particular emphasis is on comparing established
methods with a quadratic trust-region approach, which has not been widely used
for BD so far.

The remainder of the paper is structured as follows: The next section 2
generally introduces the investigated planning problem, corresponding
implementation of BD and several of its refinements. The following section 3
specifically focuses on quadratic stabilization as an additional refinement.
Subsequently, section 4 provides details on the specific planning problem and
the approach used for scenario generation that is used in section 5 to
benchmark the introduced refinements. The final section concludes and suggests
further developments.

##  2 Refined Benders algorithm

The two-stage energy planning problem decides on capacity expansion in the
first stage and operation in the second. Operation is subject to uncertainty
represented by different scenarios, for instance, corresponding to different
weather years. The following subsections first introduce the closed
formulation of the problem, followed by the basic BD implementation used to
solve it, and finally, various refinements to enhance the basic algorithm. The
notation uses lowercase letters for parameters, uppercase letters for
variables, and Greek letters for anything specific to BD.

###  2.1 Problem formulation

Eqs. 1a to 1j provide the closed formulation of the two-stage stochastic
energy planning problem. For a set of technologies $i\in I$, the problem
decides on the expansion of capacities $Exp_{y,i}$ over a set of chronological
but not necessarily consecutive years $y\in Y$. All variables of the problem
are continuous and non-negative.

According to 1d, the capacity $Capa_{y,i}$ available in each year, depends on
expansion in the current and previous years provided by $Y_{y}^{exp}$. In Eqs.
1e to 1h these capacities constrain the operation of technologies at each
time-step $t\in T$, in each scenario $s\in S$, and within each year. For
generation technologies $I^{ge}$, the capacity constraint limiting the
generation $Gen_{y,s,i,t}$ includes a scenario-dependent capacity factor
$cf_{s,i,t}$ reflecting the available share of capacity. For storage
technologies, $I^{st}$, two different capacities constrain three different
operational variables: Charged energy $St_{y,s,i,t}^{in}$ and discharged
energy $St_{y,s,i,t}^{out}$ are both constrained by the power capacity
$Capa_{y,i}^{st}$; the current storage level $St_{y,s,i,t}^{size}$ is
constrained by the energy capacity $Capa_{y,i}^{size}$. In total, net supply
from generation and storage technologies must match the exogenous demand
$dem_{y,s,t}$ in each year, time-step, and scenario, as expressed in the
energy balance in Eq. 1i. To avoid infeasibility, the balance includes a loss-
of-load variable $Lss_{y,s,t}$. The storage balance in Eq. 1d computes the
storage level at time-step $t$ based on the level at the previous time-step
$t-1$ plus charged and minus discharged energy. The storage constraint is
circular for each year, meaning the last time-step in a year is previous to
the first.

The objective function of the problem in Eq. 1a minimizes total system costs
comprised of expansion costs $U$ and the sum of operational costs $V_{y,s}$
across all years $y$ and scenarios $s$ weighted according to the scenario
probabilities $p_{s}$. Expansion costs defined in Eq. 1b depend on the
expansion variable and the specific expansion costs $c^{inv}_{y,i}$ of each
technology. Operational costs defined in Eq. 1c depend on the costs associated
with loss-of-load $c^{lss}$ and variable costs of generation technologies
$c^{var}_{y,s,i}$, for instance, fuel costs.

$\displaystyle\min$ $\displaystyle\;\;U$ $\displaystyle+$
$\displaystyle\sum_{y\in Y,s\in S}p_{s}\cdot V_{y,s}$ (1a) s.t. $\displaystyle
U$ $\displaystyle\,=\,$ $\displaystyle\sum_{y\in Y,i\in I}c^{inv}_{y,i}\cdot
Exp_{y,i}$ (1b) $\displaystyle V_{y,s}$ $\displaystyle\,=\,$
$\displaystyle\sum_{t\in T}Lss_{y,s,t}\cdot c^{lss}+\sum_{i\in
I}Gen_{y,s,i,t}\cdot c^{var}_{y,s,i}$ $\displaystyle\,\forall y\in Y,s\in S$
(1c) $\displaystyle Capa_{y,i}$ $\displaystyle\;=\;$
$\displaystyle\sum_{y^{\prime}\in Y_{y}^{exp}}Exp_{y^{\prime},i}$
$\displaystyle\;\forall y\in Y,i\in I$ (1d) $\displaystyle Gen_{y,s,i,t}$
$\displaystyle\,\leq\,$ $\displaystyle cf_{s,i,t}\cdot Capa_{y,i}^{ge}$
$\displaystyle\,\forall y\in Y,s\in S,i\in I^{ge},t\in T$ (1e) $\displaystyle
St_{y,s,i,t}^{out}$ $\displaystyle\,\leq\,$ $\displaystyle Capa_{y,i}^{st}$
$\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$ (1f) $\displaystyle
St_{y,s,i,t}^{in}$ $\displaystyle\,\leq\,$ $\displaystyle Capa_{y,i}^{st}$
$\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$ (1g) $\displaystyle
St^{size}_{y,s,i,t}$ $\displaystyle\,\leq\,$ $\displaystyle Capa_{y,i}^{size}$
$\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$ (1h) $\displaystyle
dem_{y,s,t}$ $\displaystyle\,=\,$ $\displaystyle Lss_{y,s,t}+\sum_{i\in
I^{ge}}Gen_{y,s,i,t}-\sum_{i\in I^{st}}St_{y,s,i,t}^{in}+St_{y,s,i,t}^{out}$
$\displaystyle\,\forall y\in Y,s\in S,t\in T$ (1i) $\displaystyle
St^{size}_{y,s,i,t}$ $\displaystyle\,=\,$ $\displaystyle
St^{size}_{y,s,i,t-1}+St_{y,s,i,t}^{in}-St_{y,s,i,t}^{out}$
$\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$ (1j)

The formulation in Eqs. 1a to 1j introduces the general problem structure and
all details pivotal for BD but is still stylized. The model introduced in
section 4 and used for benchmarking includes additional elements. First,
multiple regions of expansion and operation. Between regions, the model can
build transmission capacity in the first stage to exchange energy in the
second. Next, distinct energy carriers, like electricity and hydrogen, are
converted into one another by respective technologies, like hydrogen turbines
or electrolyzers. Finally, storage losses and additional constraints, imposing
upper limits on capacities or restricting operation, for instance, a yearly
emission limit. For a comprehensive formulation of the planning problem see
Göke (20212).

###  2.2 Benders decomposition

Fig. 1 is an illustrative depiction of the planning problem based on the block
structure of the constraint matrix. The expansion problem in the first stage
includes expansion and capacity variables, the constraint in 1d connecting
them and the cost definition in Eq. 1b. The second stage includes the
remaining constraints, namely the energy and storage balance, the definition
of variable costs, and the capacity restrictions. The latter links capacity
with operational variables and connects the first and second stage of the
problem.

Figure 1: Structure of closed two-stage problem

Even for a single year and scenario, the number of constraints and variables
of the expansion problem in the first stage is small compared to operation in
the second. Increasing the number of years extends both stages but increasing
the number of stochastic scenarios exclusively extends the operational
problem. Since the expansion problem is small, the total size almost linearly
depends on the number of stochastic scenarios and the problem quickly becomes
intractable, if their number increases.

For BD, the first stage expansion corresponds to the MP, while operation is
split into several independent SPs. For demonstration, consider the temporal
problem structure outlined in Fig. 2. In the first stage, two steps of
capacity expansion are modeled with a 10-year interval. Multi-year intervals
are a common modeling strategy for representing transformation pathways,
aiming to reduce computational complexity. In the second stage, there are two
operational scenarios for each expansion step. These scenarios represent
different patterns of renewable supply and energy demand. In this example,
these patterns are based on historical data for the weather years
$2000^{\prime}$ and $2008^{\prime}$.222These two years are selected at random
here for illustrative purposes. The detailed method to select representative
scenarios is introduced in section 4.2. In addition, it is conceivable to have
different operational scenarios for each expansion step, for instance to
account for changes in patterns caused by climate change.

According to the problem formulation in the previous section, operation is
modeled for each year and scenario resulting in four distinct operational
problems, as shown in Fig. 2. For a given set of capacities, each operational
problem can be solved independently, if no complicating constraint links the
variables from different SPs. An example of such a constraint is an emission
limit that is not enforced for each SP separately, but jointly across
different years and scenarios, like a carbon budget.

Figure 2: Graph of years and scenarios for exemplary problem

Fig. 3 provides the structure of the decomposed problem, analogously to Fig.
1. It shows how the operational problem is split into four independent SPs,
each consisting of different operational variables, but sharing capacities
with the MP and other SPs for the same year. Splitting the second stage
operation into several SPs is beneficial because it enables distributed
computing and prevents complexity in the second stage to scale with the number
of scenarios. In addition, it increases the number of cuts added to the MP in
each iteration improving convergence, as we will elaborate when introducing
the BD algorithm next.

Figure 3: Structure of decomposed problem

Alg. 1 presents the standard version of the BD algorithm applied in this
paper. $k$ is an iteration counter; $\zeta_{low}$ and $\zeta_{up}$ are the
current lower and upper bound of the objective function, respectively. We use
multi-cut and at each iteration add separate constraints for each sub-problem
to the MP instead of a single aggregated cut. The approach improves
convergence but increases the complexity of the MP—a favorable trade-off since
the MP is simple. To mitigate adverse effects on MP performance, the algorithm
deletes cuts that were not created or binding within the last $\eta$
iterations, with $\eta$ being a pre-set parameter. Furthermore, solving the
SPs is fully parallelized using distributed computing.

1 chose convergence tolerance $\epsilon$ and deletion threshold $\eta$

2 set $k=1$ and $\zeta_{up}\to\infty$

3 while  _$1-\frac{\zeta_{low}}{\zeta_{up}} <\epsilon$_ do

4 solve MP get $Capa_{y,i}^{(k)}$ and $U^{(k)}$

5 set $\zeta_{low}\leftarrow obj(\mathbf{MP})$

6 do in parallel

7 for  _$y\in Y,\,s\in S$_ do

8 solve $\textbf{SP}_{y,s}$ with $Capa_{y,i}^{(k)}$

9 get $V_{y,s}^{(k)}$ and $\lambda_{y,s,i}^{(k)}$

10 add $\Omega_{k,y,s}$ to MP

11 set $\delta_{k,y,s}\leftarrow k$

12

13 end for

14 for  _$l\in\\{1,\dots,k-1\\},y\in Y,\,s\in S$_ do

15 if  _$\Omega_{l,y,s}$ is binding_ then

16 set $\delta_{l,y,s}\leftarrow k$

17

18 else if  _$k-\delta_{l,y,s} >\eta$_ then

19 delete $\Omega_{l,y,s}$ from MP

20

21

22 end for

23

24 end

25

26 if  _$\zeta_{up} >U^{(k)}+\sum_{y\in Y,s\in S}p_{s}\cdot V_{y,s}^{(k)}$_
then

27 set $\zeta_{up}\leftarrow U^{(k)}+\sum_{y\in Y,s\in S}p_{s}\cdot
V_{y,s}^{(k)}$

28

29 end if

30 set $k\leftarrow k+1$

31

32 end while

Algorithm 1 Parallelized multi-cut benders, based on Conejo . (2006)

After initializing all variables, the algorithm starts iterating after line 3
solving the original MP, as described by Eqs. 2a to 2c. Since the MP only
includes the constraints of the expansion problem, its objective function
approximates the actual costs of each SP $V_{y,s}$ using the estimator
$\alpha_{y,s}$. After solving, the algorithm obtains the resulting capacities
$Capa_{y,i}^{(k)}$ and expansion costs $U^{(k)}$ for the current iteration and
sets the lower bound $\zeta_{low}$ to the current objective of the MP. In the
first iteration, $\alpha_{y,s}$ is still unconstrained and the trivial optimum
for the MP is zero meaning no capacity expansion at all.

$\displaystyle\min$ $\displaystyle\;\;U$ $\displaystyle+$
$\displaystyle\sum_{y\in Y,s\in S}p_{s}\cdot\alpha_{y,s}$ (2a) s.t.
$\displaystyle\;\;U$ $\displaystyle\,=\,$ $\displaystyle\,\,\,\,\,\sum_{y\in
Y}c^{inv}_{y,i}\cdot Exp_{y,i}$ (2b) $\displaystyle Capa_{y,i}$
$\displaystyle\;=\;$ $\displaystyle\,\sum_{y^{\prime}\in
Y_{y}^{exp}}Exp_{y^{\prime},i}$ $\displaystyle\;\forall y\in Y,i\in I$ (2c)

Next, the algorithm solves each SP in parallel with capacities fixed to the
previously computed $Capa_{y,i}^{(k)}$, as described by Eqs. 3a to 3i.
Afterward, the algorithm gets the actual costs of each SP $V_{y,s}^{(k)}$ and
the dual values $\lambda_{y,s,i}^{(k)}$ of the constraint fixing capacities in
Eq. 3i. Note that in our case, SPs are always feasible since the loss-of-load
$Lss_{y,s,t}$ serves as a slack variable for unmet demand.

$\displaystyle\min$ $\displaystyle\;V_{y,s}^{(k)}$ (3a) s.t. $\displaystyle
V_{y,s}^{(k)}$ $\displaystyle\,=\,$ $\displaystyle\sum_{t\in
T}Lss_{y,s,t}\cdot c^{lss}+\sum_{i\in I,t\in T}Gen_{y,s,i,t}\cdot
c^{var}_{y,s,i}$ $\displaystyle\,\forall y\in Y,s\in S$ (3b) $\displaystyle
Gen_{y,s,i,t}$ $\displaystyle\,\leq\,$ $\displaystyle cf_{s,i,t}\cdot
Capa_{y,i}^{ge}$ $\displaystyle\,\forall y\in Y,s\in S,i\in I^{ge},t\in T$
(3c) $\displaystyle St_{y,s,i,t}^{out}$ $\displaystyle\,\leq\,$ $\displaystyle
Capa_{y,i}^{st}$ $\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$
(3d) $\displaystyle St_{y,s,i,t}^{in}$ $\displaystyle\,\leq\,$ $\displaystyle
Capa_{y,i}^{st}$ $\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$
(3e) $\displaystyle St^{size}_{y,s,i,t}$ $\displaystyle\,\leq\,$
$\displaystyle Capa_{y,i}^{size}$ $\displaystyle\,\forall y\in Y,s\in S,i\in
I^{st},t\in T$ (3f) $\displaystyle dem_{y,s,t}$ $\displaystyle\,=\,$
$\displaystyle Lss_{y,s,t}+\sum_{i\in I^{ge}}Gen_{y,s,i,t}-\sum_{i\in
I^{st}}St_{y,s,i,t}^{in}+St_{y,s,i,t}^{out}$ $\displaystyle\,\forall y\in
Y,s\in S,t\in T$ (3g) $\displaystyle St^{size}_{y,s,i,t}$ $\displaystyle\,=\,$
$\displaystyle St^{size}_{y,s,i,t-1}+St_{y,s,i,t}^{in}-St_{y,s,i,t}^{out}$
$\displaystyle\,\forall y\in Y,s\in S,i\in I^{st},t\in T$ (3h) $\displaystyle
Capa_{y,i}$ $\displaystyle=$ $\displaystyle
Capa^{(k)}_{y,i}\,:\,\lambda_{y,s,i}^{(k)}$ $\displaystyle\forall i\in I$ (3i)

Next, the Benders cuts $\Omega_{k,y,s}$ defined by Eq. 4 are computed and
added to the MP. Each cut improves the estimator $\alpha_{y,s}$ based on the
the exact result of the SP $V_{y,s}^{(k)}$ and its gradient
$\lambda_{y,s,i}^{(k)}$ at the point $Capa_{y,i}^{(k)}$. In other words, by
adding hyperplanes restricting the estimator, BD performs a piece-wise linear
approximation of the unknown but convex function that assigns capacities to
SP-costs (Conejo ., 2006). In the literature, cuts of this form are termed
optimality cuts, as opposed to feasibility cuts that are added when the SP was
infeasible, which cannot occur in our case.

$\Omega_{k,y,s}:\alpha_{y,s}\geq V^{(k)}_{y,s}+\sum_{i\in
I}\lambda_{y,s,i}^{(k)}\cdot(Capa_{y,i}-Capa^{(k)}_{y,i})$ (4)

Afterward, the value $\delta_{k,y,s}$ tracking in which iteration a cut was
created or binding the last time is initialized based on the number of the
current iteration $k$. Parallel to solving the SPs, lines 13 to 18 of the
algorithm loop over all previous cuts updating $\delta_{k,y,s}$ and deleting
cuts that were not created or binding within the last $\eta$ iterations.

Once all SPs are solved, the total costs at $Capa_{y,i}^{(k)}$ are computed by
summing expansion costs $U$ and accurate operational costs $V_{y,s}^{(k)}$. If
these undercut the current $\zeta_{up}$, $Capa_{y,i}^{(k)}$ is a new best
solution and $\zeta_{up}$ is adjusted accordingly. Finally, the iteration
counter increases by one and the next iteration starts, if the optimality gap
does not exceed the predefined tolerance yet.

###  2.3 Algorithmic refinements

Since the basic BD algorithm converges slowly, there are various methods to
improve its performance. Most refinements however aim at problems with a
complex combinatorial MP, rather than large SPs, like the two-stage planning
problem examined in this paper.

Several refinements have been considered and to some extent tested, but their
full implementation was ultimately discarded for various reasons. We did not
deploy advanced cut generation, because it requires solving an adapted version
of the dual SP in addition to the primal SP, which is not beneficial if the
SPs are large. Another approach is to utilize faster surrogate models instead
of the actual SP. For instance, Mazzi . (2021) employ a fitted saddle function
that is trained to approximate the SP from above and thereby generates valid
Benders cut. However, we have observed that such methods only achieve
sufficient accuracy when the problem has a large number of small SPs and the
dimensionality of in- and outputs for the surrogate model is low, implying
that the number of complicating capacity variables is small. Other
implementations adopt a strategy of reducing the number of scenarios
considered in the initial iterations. However, this approach is primarily
applicable when dealing with a large number of scenarios that cannot be
completely parallelized and are randomly sampled from a stochastic model. In
our specific case, as described in the section (Section 4), the scenarios are
based on empirical data from different historical years resulting in a
comparatively small number of SPs that can be fully parallelized. Finally, we
did not apply any further decomposition, for example, Dantzig–Wolfe to address
the complicating storage constraints in Eq. 3h of the SPs, or dual dynamic
programming to solve the multi-year expansion problem sequentially (Sepulveda,
2020). The existing research indicates that additional decomposition of the
already parallelized LP will not result in any improvement in computation
time.

Several refinements to the original algorithm have already been introduced in
the previous subsection 2.2. These include the multi-cut reformulation,
parallelization of SPs, and the procedure to delete non-binding cuts after a
predefined number of iterations. In the literature, strategies like the latter
to manage the bundle of subgradients, or cuts, are also known as
$\beta$-strategies.333This $\beta$ is not to be confused with the parameter of
the level bundle method discussed in 3 They generally reduce the solution time
of the MP but can impair convergence of the overall algorithm (Frangioni,
2002).444See the beginning of section 3 for context on bundle methods.

Furthermore, the outlined algorithm is subject to several practical
refinements. All capacity variables have an upper limit that is far beyond any
reasonable value to pre-limit the solution space of the MP. To improve
numerical properties, all equations are scaled according to the methodology
outlined in (Göke, 20211). If scaling does not suffice, small values of
$Capa_{y,i}^{(k)}$ in the SPs are rounded off to zero. The same applies to
$\lambda_{y,s,i}^{(k)}$ when adding cuts to the MP. For storage, we enforce an
upper and lower bound on the ratio between storage and energy capacity that
reflects technical restrictions and confines the solution space.

The following subsections provide a detailed description of two implemented
refinements: inexact cuts and valid inequalities. The subsequent section
discusses stabilization methods separately.

####  2.3.1 Inexact cuts

According to Zakeri . (2000), Benders cuts in Eq. 4 are not required to be
derived from an optimal solution and can be based on any feasible solution to
the SP. Although inexact cuts can have an adverse effect on convergence, they
have the—potentially compensatory—benefit of reducing the solution time of the
SPs.

In our implementation of the BD algorithm, we utilize the parameters of the
Barrier algorithm solving the individual SPs to compute inexact cuts. Being an
interior point solver, early termination of barrier will provide a non-optimal
solution that is within the interior of the solution space and thereby
feasible. Especially in early iterations, when the MP solution is still far
from the optimum, inexact cuts are favorable and have less impact on
convergence. Accordingly, the convergence tolerance of the SPs solved with
Barrier decreases with the optimality gap of the overall algorithm. Fig. 4
presents three different interpolation methods for this process. In all cases,
the algorithm starts at a tolerance of $1e-2$ and decreases to $1e-8$, but the
reduction either follows an exponential, linear, or logarithmic curve. As a
result, each curve imposes a distinct trade-off between the convergence of the
overall algorithm and the solution time of the SPs.

Figure 4: Interpolation methods for convergence tolerance of SPs

In all cases, solving the SPs does not include the crossover phase of Barrier
that eventually determines an optimal basis solution. Crossover typically has
a negligible impact on results, but significantly affects computation time.

####  2.3.2 Valid inequalities

Valid inequalities are constraints added to the MP that are redundant in the
closed formulation of the problem but effectively reduce the solution space of
the MP (Cordeau ., 2006). Thus, valid inequalities aim to enhance the
convergence of BD at the expense of increasing the MP and its solution time.
Often valid inequalities are derived from the SPs and avoid their
infeasibility.

In examined two-stage problem, the MP does not include an energy balance and
will frequently underinvest in generation capacity as a result, in particular
in early iterations. To address this, Wang . (2013) add operational
constraints, including the energy balance, of one scenario to the MP. In our
case, this approach is not viable due to the size of the individual scenarios
and corresponding SPs. Instead, we add aggregated operational constraints in
Eqs. 5a and 5b to the MP for each scenario.

$\displaystyle\sum_{t\in T}dem_{y,s,t}$ $\displaystyle\,\leq\,$
$\displaystyle\sum_{i\in I^{ge}}Gen_{y,s}$ $\displaystyle\,\,\,\,\forall y\in
Y,s\in S$ (5a) $\displaystyle Gen_{y,s,i}$ $\displaystyle\,\leq\,$
$\displaystyle\sum_{t\in T}cf_{s,i,t}\cdot Capa_{y,i}^{ge}$
$\displaystyle\,\,\,\,\forall y\in Y,s\in S,i\in I^{ge}$ (5b)

The first equation 5a enforces a yearly energy balance setting a lower bound
for generation that is achieved when hourly generation patterns precisely
match demand, without relying on storage methods that incur energy losses. The
second equation 5b establishes a connection between the lower limit on
generation and the capacity variables. As discussed in section 2.1, the
problem formulation in the paper is still stylized and omits the multiple
regions in the case study connected by transmission infrastructure. Therefore,
the applied valid inequalities are additionally aggregating by region assuming
lossless and unrestricted transmission as a lower bound.

As highlighted in section 2.1, the problem formulation in the paper remains
stylized and does not account for the multiple regions connected by
transmission infrastructure in the case study. Consequently, the applied valid
inequalities are also aggregating demand and generation by region, assuming
lossless and unrestricted transmission as a lower bound.

##  3 Quadratic stabilization

Quadratic stabilization is the pivotal refinement of the developed BD
algorithm. The stabilization methods discussed below are not restricted to
stochastic programming but derive from a much broader field of literature
concerned with cutting plane methods. Cutting plane methods iteratively
minimize convex functions, which may be non-differentiable, by approximating
the function through outer linearizations (Kelley, 1960). Consequently, BD is
a cutting plane method specifically for two-stage problems.

We begin by representing the BD algorithm above in a more general form to
connect it to the literature on stabilization in non-smooth convex
optimization.

Let $\phi_{y,s}(Capa)$ be the optimal objective value of the SP for investment
year $y$ and scenario $s$ as a function of complicating variable choice
$Capa$. This function as well as the weighted average of such functions can be
shown to be polyhedral (and thus convex and non-smooth) (Shapiro ., 2009,
Prop. 2.2 & 2.3), which allows constructing a cutting plane model.

The cutting plane approximation of $\phi_{y,s}(Capa)$ at iteration $k$
corresponds to the minimum of the estimator $\alpha_{y,s}$ restricted by the
cuts $\Omega_{k^{\prime},y,s}$, as shown in Eq. 6.555This is equivalent to the
piece-wise maximum of all linearizations across the domain of $Capa$, as the
cutting plane model is usually defined in the convex optimization literature.

$\displaystyle\hat{\phi}_{y,s}^{(k)}(Capa)=\displaystyle{\min_{k>k^{\prime}}\\{\alpha_{y,s}\,|\,\Omega_{k^{\prime},y,s}\\}}$
(6)

The sum of approximations over all SPs $\hat{\phi}^{(k)}(Capa)$ can
accordingly be written as Eq. 7.

$\displaystyle\hat{\phi}^{(k)}(Capa)=\sum_{y\in Y,s\in
S}p_{s}\hat{\phi}_{y,s}^{(k)}(Capa)$ (7)

By sending a candidate solution $Capa^{(k)}$ to an oracle, which corresponds
to solving the SPs in our case, we obtain
$\sum_{y,s}p_{s}\phi_{y,s}(Capa^{(k)})$ and subgradients
$\lambda_{y,s,i}^{(k)}$, which are in the subdifferentials of the respective
$\phi_{y,s}(Capa^{(k)})$.666Subgradients generalize the concept of gradients
to non-differentiable functions. $g$ is a subgradient of $f(x)$ if $f(y)\geq
f(x)+\langle g,y-x\rangle\ \forall y$. The subdifferential of $f(x)$ is the
set of all subgradients at $x$. If the function is differentiable at $x$, the
subdifferential is a singleton. On this basis, we improve the cutting plane
model in Eq. 7 by adding new cuts and then solve the MP to obtain a new
candidate solution $Capa^{(k+1)}$.

Although the proof of convergence for $\min\hat{\phi}(Capa)\to\min\phi(Capa)$
is simple (Kelley, 1960), cutting plane algorithms suffer from two
computational issues. First, the algorithm is inherently unstable since the
capacities $Capa$ in two consecutive iterates can be arbitrarily far apart
extreme points of the feasible space. In our case, this problem is
particularly pronounced, because the feasible region is hardly restricted and
completely continuous. In fact, even enforcing the optimum as a starting value
cannot be expected to improve the solution time, since the algorithm is
initially unable to evaluate the quality of the solution. Second, the size of
the model approximating $\phi(Capa)$, in the case of BD the number of cuts,
can become prohibitively large as the number of iterations $k$ increases and
slow the algorithm (Frangioni, 2020). It should be emphasized that these
issues are mutually reinforcing: Oscillation increases the number of
iterations, which in turn slows down the convergence of the approximating
model $\hat{\phi}^{(k)}(Capa)$.

The term ”bundle methods” encompasses methods to address these two issues by
managing the ”bundle” of information accumulated throughout the iterations
(Frangioni, 2002). In the following, our exclusive focus will be on
stabilization methods that center each iteration around a reference point in
order to prevent oscillation (Bonnans ., 2006). Since these methods improve
convergence at the expense of adding complexity to the MP, they are well-
suited to the examined problem with a simple MP but costly oracles, viz. large
individual SPs. The second issue of an increasingly complex approximation is
already addressed by the previously described method to remove non-binding
cuts after a certain number of iterations. Additionally, stabilization reduces
the overall number of iterations, thereby mitigating problems associated with
excessive growth of the bundle as well.

The literature on bundle methods also considers the combination of
stabilization and inexact oracles (van Ackooij ., 2016; van Ackooij de
Oliveira, 2014; de Oliveira ., 2014; de Oliveira Sagastizábal, 2014), in which
$\lambda^{(k)}_{y,s}$ and/or $\phi(Capa^{(k)})$ are not determined exactly. In
fact, the use of inexact cuts described in section 2.3.1 is an example of
asymptotically exact oracles, which can readily be combined with the bundle
methods below (Frangioni, 2020, Ch. 3.5). Other inexact methods would require
further modification of the applied stabilization methods.

The following section introduces the common characteristics of stabilization
methods. Next, we describe the proximal bundle and level bundle methods, the
two most common techniques for quadratic stabilization in the literature, and
explain how they have been implemented in the developed BD algorithm.
Afterward, we present a theoretically equivalent but not widely deployed
trust-region approach as an alternative in greater detail.

###  3.1 Characteristics of stabilization methods

The convex optimization literature has proposed three main methods to address
the instability of cutting plane methods that share certain characteristics.

These stabilization methods revolve around a stability center, a candidate
solution to the MP, and aim to keep the next iterate solution relatively close
to the current center. In theory, various measures are conceivable to quantify
the distance between the stability center and the next iterate, but the
corresponding literature typically relies on the Euclidean norm. At least for
the proximal bundle method, this has been justified by the effective
construction of a ”poorman’s Hessian”, which approximates valuable second-
order information otherwise absent in the polyhedral setting (Frangioni,
2020). In addition, modern solvers can handle optimization problems with
quadratic equations comparatively well. The stability center also enables to
benefit from initial heuristic solutions when using stabilization methods for
cutting plane methods.

Typically, the stability center is defined in terms of the input to the
optimal value function, in our case the capacity vector $Capa$. However, the
stability center can be formulated in terms of any MP variables that are
implicitly defined by the input of the original value function, although
selecting variables closely connected to the input of the optimal value
function appears sensible. In our case, we opt for selecting the expansion
variable instead, as many capacities are fixed in practical applications.
Consequently, including them in the center would introduce unnecessary
complexity to the stabilization method.

An update of the stability center is called a serious step and an iteration
that keeps the current stability center is called a null step. Contrary to
what the name suggests, a null step is still beneficial since it improves the
cutting plane model locally. The literature offers various methods to
distinguish a null step and a serious step. For example, a serious step can be
defined as a step improving the objective by a specified value compared to the
current stability center. In an alternative definition, the ratio between the
actual improvement determined by $\phi(Capa)$ and the improvement estimated by
$\hat{\phi}(Capa)$ must exceed a certain threshold. In our case, stabilization
achieved the best results when the algorithm makes a serious step in each
iteration improving the current best solution, even if just marginally.

To control the distance between the current stability center and the next
iterate, all stabilization methods deploy a dynamic parameter. Although
convergence is even guaranteed for a static parameter, updating the parameter,
also referred to as the t-strategy, can substantially improve performance
(Frangioni, 2002). It is important to point out that all the discussed
stabilization methods are theoretically equivalent. This means for a given
stability center, there is a parameterization for each method that will give
the same candidate for $Capa$ in the next iteration (Frangioni, 2002).
However, determining this parameter for each method is difficult, and in
practice, the t-strategy greatly impacts the performance of the algorithm (de
Oliveira Solodov, 2016). Stabilization with a poor t-strategy can even be
inferior to an unstabilized cutting plane algorithm. For instance, an
overstabilized MP, will not benefit from the global nature of the cutting
plane model and potentially results in more iterations than an unstabilized
cutting plane method (cf. Frangioni, 2020, Fig. 3.2).

###  3.2 Proximal bundle and level bundle methods

The proximal bundle and level bundle method presented in this section are the
two most common quadratic stabilization methods for cutting plane algorithms.
For the comparison in this paper, we implement effective but standard versions
of these algorithms, but it should be noted that the literature also presents
further variations and more elaborate t-strategies (Rey Sagastizábal, 2002;
Schramm Zowe, 1992; de Oliveira Solodov, 2016).

Proximal bundle methods add a regularizing term to the MP penalizing
deviations from the current stability center (Kiwiel, 1990; Frangioni, 2002).
Let the stability center based on previous iterations be $Exp^{up}$, then the
next iterate $Capa^{(k)}$ is computed according to Eq. 8, an adapted version
of the MP in Eqs. 2a to 2c. In the equation, $U(Capa)$ refers to the expansion
costs according to Eq. 2b.

$\displaystyle
Capa^{(k)}=\arg\min\left\\{U(Capa)+\hat{\phi}^{(k)}(Capa)+\mu^{(k)}\|Exp-
Exp^{up}\|_{2}^{2}:Capa_{y,i}=\sum_{y^{\prime}\in
Y_{y}}Exp_{y^{\prime},i}\right\\}$ (8)

The dynamic parameter $\mu^{(k)}$ penalizes deviations from the current
stability center. Determining an effective t-strategy for this parameter is
not trivial (de Oliveira Solodov, 2016). Early proposals of the algorithm use
a constant term but most implementations double the penalty after a null step
and half the penalty after a serious step (Ruszczyński, 1986; Zverovich .,
2012). In addition, an upper bound is typically imposed on the penalty. Our
implementation follows this approach. In addition, it resets the penalty to
its starting value if it was at its upper bound for two consecutive
iterations, since we observed poor convergence when the penalty is at the
upper bound for multiple iterations.

Level bundle methods do not add a penalty to the original objective function,
but directly redefine the objective to minimize the distance between the
stability center and the next iterate (Lemaréchal ., 1981). An additional
constraint sets an upper limit $\ell^{(k)}$, the dynamic parameter of this
method, for the estimated objective value forcing the problem to move away
from the center (Lemaréchal ., 1995), as shown in Eq. 9.

$\displaystyle Capa^{(k)}=\arg\min\left\\{\|Exp-
Exp^{up}\|_{2}:U(Capa)+\hat{\phi}^{(k)}(Capa)\leq\ell^{(k)},\
Capa_{y,i}=\sum_{y^{\prime}\in Y_{y}}Exp_{y^{\prime},i}\right\\}$ (9)

The typical choice for $\ell^{(k)}$ is the linear combination in Eq. 10
controlled by the parameter $\beta$ between the objective at the current
stability center, an upper limit for the problem, and the current objective of
the unstabilized MP, a lower limit for the problem.

$\displaystyle\ell^{(k)}=$
$\displaystyle\beta\cdot(U(Capa^{up})+\phi^{(k)}(Capa^{up}))$ (10)
$\displaystyle+(1-\beta)\cdot\min\left\\{U(Capa)+\hat{\phi}^{(k)}(Capa):Capa_{y,i}=\sum_{y^{\prime}\in
Y_{y}}Exp_{y^{\prime},i}\right\\}$

A common problem with level bundle methods is that the stabilized MP can be
infeasible if the imposed level constraint is overly ambitious. In our
implementation, we resolve this problem by replacing the lower value in Eq. 10
with the current level parameter $\ell^{(k)}$ until the problem is feasible.

The benefit of defining a hyperparameter $\ell^{(k)}$, by means of $\beta$, is
that its value is directly linked to the objective function. Setting the
starting value for the dynamic parameter of proximal bundle methods is—at
least in the standard version described here—rather an arbitrary decision.

###  3.3 Trust-region method

Trust-region methods leave the objective function of the MP unchanged but
restrict its solution space instead.

For continuous problems, most applications consider the seminal version with
an $l_{\infty}$ norm in order to restrict the distance between the stability
center and the next iterate, resulting in a linear program (Marsten ., 1975).
We tested this method but it was, in line with the existing literature (e.g.
Frangioni Gorgone, 2014), not found to stabilize the algorithm effectively.
For binary variables, trust-regions can use the $l_{0}$-norm, or Hamming
distance, and remain linear Santoso . (2005).

Again, we only consider the Euclidean norm as a distance measure resulting in
a quadratic trust-region. Accordingly, the proximal bundle method from the
previous section can be considered the Lagrangian relaxation of the quadratic
trust-region method (Frangioni, 2020, Ch. 3.3.2), highlighting the theoretical
equivalence. Recent improvements of commercial solvers made the addition of a
quadratic constraint to a linear problem more viable. To the best of our
knowledge, quadratic trust regions for continuous problems have been mentioned
but not applied in the existing literature.

To enforce the trust-region, Eqs. 11 and 12 are added to the MP in in Eqs. 2a
to 2c. To link the radius with the decision variables, it is defined in
relation to the $l^{1}$-norm of the current stability center $Exp^{up}$. If
the algorithm makes a serious step and the stability center is updated, the
radius is adjusted accordingly.

$\displaystyle\;||Exp-Exp^{up}||_{2}$ $\displaystyle\;\;\leq$
$\displaystyle\psi$ $\displaystyle\cdot$ $\displaystyle||Exp^{up}||_{1}$ (11)
$\displaystyle\iff$ $\displaystyle\sum_{y\in Y,i\in I}$
$\displaystyle(\,Exp_{y,i}-Exp_{y,i}^{up}\,)^{2}$ $\displaystyle\;\;\leq$
$\displaystyle\;\;(\,$ $\displaystyle\psi$ $\displaystyle\cdot$
$\displaystyle\sum_{y\in Y,i\in I}Exp_{y,i}^{up}\,)^{2}$
$\displaystyle\;\;\;\;\forall y\in Y,s\in S$ (12)

The parameter $\psi$ that controls the ratio between the $l^{1}$-norm of the
current stability center is the dynamic parameter of the trust-region method.
If the trust-region is not a binding constraint for the MP, the relative
radius $\psi$ is halved. Since the problem is linear, the trust-region is only
non-binding if the optimum is already within it, and reducing the radius is
sensible.

In the following, we illustrate the trust-region stabilization based on four
consecutive iterations for a stylized example problem in Fig. 5. It includes
three unbounded expansion variables, restricted by a spherical trust-region
according to Eq. 12. The stability center is indicated in orange,
corresponding to the current best, and the solution at the current iteration
$Exp^{(k)}$ in yellow or red, depending on whether it improves the current
best or not. The minimizer of $\phi(Capa)$ is shown in green. The first
iteration in (a) results in an objective above the current best, but the
trust-region is binding and remains unchanged. Due to the added cuts, the
second iteration (b) improves the current best and the trust-region is re-
centered for the third iteration. Therefore, the radius is also adapted based
on the $l^{1}$-norm of the new solution. The third iteration in (c) does not
improve the current best, but the trust-region is not binding anymore, and its
radius halved accordingly. The fourth iteration with a smaller trust-region in
(d) does not improve the solution.

Figure 5: Adaption of trust-region during iteartion

Alg. 2 extends the BD algorithm introduced in the previous section according
to the dynamic trust-region approach. The adjustments for other stabilization
methods are analogous and can be found in the cited literature.

If the algorithm uses initialization, it first solves the closed problem (CP)
for the scenario with the highest probability $p_{s}$. Afterward, all SPs are
solved in parallel for $Capa_{y,i}^{(0)}$ to set $\zeta_{up}$, and the
corresponding Benders cuts are added to the MP. Without initialization, the
steps prior to iteration do not change.

During iteration, the algorithm now solves the quadratic master-problem (QMP)
including the trust-region to determine $Capa_{y,i}^{(k)}$. Parallel to
solving the SPs for $Capa_{y,i}^{(k)}$, the regular MP is solved as well to
determine its objective value and expansion costs. The MP’s objective then
sets the lower bound $\zeta_{low}$ since the trust-region distorts the QMP’s
objective. The expansion costs of the MP and QMP, $U^{(k)}_{std}$ and
$U^{(k)}_{qua}$ respectively, determine if the trust-region is binding and the
relative trust-region radius $\psi$ is halved in line 42. In terms of coding,
this approach is more convenient than explicitly computing the slack of the
actual constraint. The trust-region is only halved, if the current iteration
did not change the current best solution $Exp_{y,i}^{up}$. If the iteration
did improve the current best, the stability center is updated in line 39 and
the radius of the trust-region is adapted. As a result, the absolute radius of
the trust-region can decrease making the updated trust-region binding again.

1 chose convergence tolerance $\epsilon$ and deletion threshold $\eta$

2 chose initial relative radius $\psi$

3 if  _initialize trust-region $\mathit{true}$_ then

4 solve CP with $S^{\prime}\coloneqq\\{s\in S\,|\,p_{s}=\max{(p_{s})}\\}$ and
get $Exp_{y,i}^{(0)}$

5 do in parallel

6 for  _$y\in Y,\,s\in S$_ do

7 solve $\textbf{SP}_{y,s}$ with $Capa_{y,i}^{(0)}$

8 get $V_{y,s}^{(0)}$ and $\lambda_{y,s,i}^{(0)}$

9 add $\Omega_{0,y,s}$ to QMP and MP

10 set $\delta_{0,y,s}\leftarrow 0$

11

12 end for

13

14 end

15

16 set $\zeta_{up}\leftarrow U^{(0)}+\sum_{y\in Y,s\in S}p_{s}\cdot
V_{y,s}^{(0)}$

17 set $Exp_{y,i}^{up}\leftarrow Exp_{y,i}^{(0)}$ and update trust-region of
QMP

18

19else

20 set $\zeta_{up}\to\infty$

21 end if

22set $k=1$

23 while  _$1-\frac{\zeta_{low}}{\zeta_{up}} <\epsilon$_ do

24 solve QMP and get $Capa_{y,i}^{(k)}$, $Exp_{y,i}^{(k)}$ and $U^{(k)}_{qua}$

25 do in parallel

26 for  _$y\in Y,\,s\in S$_ do

27 solve $\textbf{SP}_{y,s}$ with $Capa_{y,i}^{(k)}$

28 get $V_{y,s}^{(k)}$ and $\lambda_{y,s,i}^{(k)}$

29 add $\Omega_{k,y,s}$ to QMP and MP

30 set $\delta_{k,y,s}\leftarrow k$

31

32 end for

33 for  _$l\in\\{1,\dots,k-1\\},y\in Y,\,s\in S$_ do

34 if  _$\Omega_{l,y,s}$ is binding_ then

35 set $\delta_{l,y,s}\leftarrow k$

36

37 else if  _$k-\delta_{l,y,s} >\eta$_ then

38 delete $\Omega_{l,y,s}$ from MP

39

40

41 end for

42 solve MP and get $U^{(k)}_{std}$

43 set $\zeta_{low}\leftarrow obj(\mathbf{MP})$

44

45 end

46

47 if  _$\zeta_{up} >U^{(k)}+\sum_{y\in Y,s\in S}p_{s}\cdot V_{y,s}^{(k)}$_
then

48 set $\zeta_{up}\leftarrow U^{(k)}+\sum_{y\in Y,s\in S}p_{s}\cdot
V_{y,s}^{(k)}$

49 set $Exp_{y,i}^{up}\leftarrow Exp_{y,i}^{(k)}$ and update trust-region of
QMP

50

51 else

52 if  _$U^{(k)}_{std}=U^{(k)}_{qua}$_ then

53 set $\psi\leftarrow 0.5\cdot\psi$ and update trust-region of QMP

54

55 end if

56

57 end if

58 set $k\leftarrow k+1$

59

60 end while

Algorithm 2 Parallelized multi-cut benders with dynamic quadratic trust-region

##  4 Case study

This section introduces the case-study used to test and benchmark the
algorithms introduced previously. The ambition of the case-study is not to
achieve highly accurate system planning, but to provide a plausible problem
covering all critical elements of renewable energy systems. All files to re-
produce the case-study are placed in a public repository linked in the
Supplementary material section.

###  4.1 Applied planning model

The applied planning model includes two expansion years, 2030 and 2040, and
focuses on the power sector. To cover a diverse range of system setups, the
planned system is not subject to any emission constraints in 2030 but must
fully decarbonize by 2040. Spatially the model covers four distinct regions:
France, Belgium, the Netherlands, and Germany.

Figure 6 introduces all considered technologies, depicted as gray circles, and
their interaction with energy carriers, depicted as colored squares. Exogenous
demand in the model is limited to the carrier electricity which is modeled at
an hourly resolution for the entire year. Hydrogen uses a daily resolution;
oil and gas a yearly resolution.

Figure 6: Graph of model elements

In the graph, entering edges of technologies refer to their input carriers;
outgoing edges relate to outputs. For example, the generation technology
electrolysis uses electricity as an input to generate hydrogen. Storage
technologies, like pumped storage, have an entering and an outgoing edge to
represent charging and discharging. Beyond pumped storage, the model includes
batteries for short-term and hydrogen tanks for long-term energy storage.
Including these technologies is critical due to their importance for
decarbonized systems and their impact on the structure of the planning
problem.

Capacities for wind and PV are subject to an upper limit that restricts
expansion. Pumped storage, run-of-river, and hydro reservoirs have capacities
fixed to today’s levels. Furthermore, hydro reservoirs are modeled as storage
with an exogenous inflow instead of flexible charging. Beyond the listed
technologies, the model also decides on the expansion and operation of
transmission to exchange electricity and hydrogen between regions.

###  4.2 Scenarios for weather years

To compare the introduced variations of BD, we solve the model as a two-stage
planning problem with multiple scenarios. Each scenario represents a weather
year with different patterns and total levels of electricity demand and
capacity factors for wind and PV. For this purpose, we build on corresponding
data for 39 weather years from 1980 to 2018 published in H. Bloomfield .
(2020), based on re-analysis data from the NASA’s MERRA-2 dataset.

Stochastic programming often reduces $n$ scenarios of historical data, in our
case 39 weather years, to a smaller number $m<n$ that captures the empirical
distribution but improves computational tractability (Dupačová, J. and Gröwe-
Kuska, N. and Römisch, W., 2003; Rujeerapaiboon ., 2022). Scenario reduction
methods determine a set of representative scenarios $S$ and a set of
corresponding weights $\\{p_{s}:s\in S,\sum_{s\in S}p_{s}=1,p_{s}\geq 0\\}$,
such that the weighted scenario average of second stage costs approximates the
sample average over all $n$ scenarios.777The sample average, in turn, is
assumed to converge asymptotically to the expected cost of the second stage,
provided a suitable law of large numbers holds. There is a variety of
reduction methods for this purpose ranging from naïve sampling to moment
matching or clustering methods (Kaut, 2021). The choice of scenario reduction
method and the optimal number of scenarios $m$ for adequate planning of
renewable systems is an open question that is beyond the scope of this paper.
For a meaningful comparison of solution methods, as conducted in Section 5
below, it is sufficient to ensure the $m$ selected scenarios are not extremely
similar and bias the algorithmic performance as a result.

In this paper, we used a problem-dependent k-mediod algorithm, as proposed by
Bertsimas Mundru (2022), for selecting $m$ representative scenarios and
corresponding probability weights. The method is advantageous for large-scale
planning problems since it reduces the high-dimensional clustering problem
including hourly data on demand and capacity factors across multiple countries
to a much smaller clustering problem in $\mathbb{R}_{+}$, the space of system
costs. To perform the clustering, we calculate a matrix of symmetric distances
between scenarios, or weather years, where the $s$-th row and $s^{*}$-th
column is computed according to Eq. 13:

$\displaystyle
d_{s,s^{*}}=\frac{1}{2}\left[Z_{s,s^{*}}-Z_{s,s}+Z_{s^{*},s}-Z_{s^{*},s^{*}}\right]$
(13)

where $Z_{s,s^{*}}$ refers to total system costs when the deterministic
problem for scenario $s$ is solved with capacities fixed to results of the
deterministic problem for scenario $s^{*}$. For instance, $Z_{1980,2008}$
corresponds to total system costs, including infeasibility costs, when solving
the deterministic model for the weather year 1980 with capacities originally
computed for 2008. The clustering algorithm selects $m$ representative
scenarios and assigns each of the $n-m$ scenarios remaining in the sample to
one of the representative scenarios, such that the sum of distances between
representative scenarios and assigned scenarios as defined above, is
minimized.

Fig. 7 illustrates the results of the clustering algorithm for $m=6$. Each
color represents a group and the node of the representative scenario for this
group is equipped with its weight $p_{s}$, which reflects the relative group
size. The color and thickness of connecting lines between the nodes indicate
similarity, the inverse of the distance $d_{s,s^{*}}$ as defined above.

Figure 7: Similarity graph for problem-dependent scenario reduction with $m=6$
clusters

##  5 Results

This section first compares the different stabilization methods and additional
refinements. Subsequently, the algorithm’s optimal configuration is
benchmarked by comparing its performance to an off-the-shelf solver for the
deterministic equivalent. The comparison is based on solving the case-study
while varying the number of stochastic scenarios from 2 to 16. The MP and each
SP run on independent nodes of a computing cluster, each with four cores and
16 GB of memory.888Core models can differ between runs and nodes. To ensure
hardware differences do not bias results, conspicuous performance outliers
were re-run always reproducing the initial results. Accordingly, the number of
nodes amounts to one for the MP plus twice the scenario number for the SPs,
since the case-study covers two expansion years, 2030 and 2040.

For solving problems, each node deploys the Barrier implementation of Gurobi
without Crossover and the NumericFocus parameter set to zero. For BD, we
delete unused cuts after 20 iterations and choose a convergence tolerance
$\epsilon$ of 1%. At smaller tolerances, the algorithm sometimes fails to
converge due to rounding (see the end of section 2.2). All deployed variations
of BD are implemented in the version of the AnyMOD.jl modeling framework
linked in the Supplementary material.

###  5.1 Comparison of stabilization methods

First, we compare the three stabilization methods from section 3 for the case
with two and four scenarios. For the proximal bundle method, we vary the
starting value of the dynamic parameter $\mu_{(1)}$ between 0.05 and 0.5. In
addition, it is subject to an upper limit of 5. An upper limit of 50 was
tested as well but did not achieve better results. For the level bundle
method, we vary the $\beta$ parameter between $0.125$ and $0.5$. Finally, we
vary the relative starting radius $\psi_{(1)}$ for the trust-region method
between 0.005 and 0.02. In all cases, the stability center is initialized with
the solution of the deterministic problem for the scenario with the highest
probability. Furthermore, all runs for the comparison use inexact cuts with
the logarithmic interpolation method, but no valid inequalities (see section
2.3). The same applies to results presented in this section, if not stated
otherwise.

The results of the comparison in Fig. 8 demonstrate the substantial speed-up
the trust-region approach achieves over the standard algorithm decreasing the
solution time by a factor of four to ten. The acceleration results from
greatly reducing the number of iterations, which outweighs the increase in
time per iteration due to the added quadratic equation. The results show that
the performance of all stabilization methods is roughly comparable. However,
in the case of the proximal bundle method and two scenarios, convergence slows
significantly in later iterations and eventually this variant is slower than
the unstabilized algorithm. While the proximal bundle method exhibits faster
convergence with four scenarios, it does not demonstrate any advantages over
the other methods in this case either. Therefore, it was excluded from further
consideration. For the level bundle method, further analysis will use a
$\beta$ parameter of 0.25 since it yielded the best results; for the trust-
region method a relative starting radius $\phi_{(1)}$ is used for the same
reason.

Figure 8: Comparison of stabilization methods and no stabilization (black
dots)

Fig. 9 shows the impact the initialization has on the performance of the
stabilized algorithm. Due to the time spent on computing the heuristic
solution, with initialization, the first iteration starts later, and
accordingly, points are offset to the right. However, in all cases, the
initialization significantly compensates for this delay and greatly
accelerates the convergence of the algorithm.

Figure 9: Impact of initialization on stabilization method

Furthermore, results indicated that adding storage technologies has a much
greater effect on iterations and solution time than adding generation
technologies. In contrast to generation, storage is not constrained by a
single capacity variable, but multiple interlinked capacities and their
interaction are difficult to capture for BD. For instance, if only the
capacity constraints on storage output are binding in a SP, the shadow
variable $\lambda_{y,s,i}^{(k)}$ is large for power output, but zero for
energy capacity. Consequently, the algorithm will solely expand output
capacity in the following iteration and not consider how the energy capacity
constraints operation as well.

###  5.2 Impact of other refinements

In the next step, we analyze the impact other refinements have on the
algorithm’s performance.

First, Table 1 shows results for the different interpolation methods
controlling the convergence tolerance of the SPs in dependence of the MP’s
optimality gap, as introduced in section 2.3.1. The results indicate that
limiting small tolerances to later iterations with a small optimality gap is a
sensible strategy to reduce computation time. With the logarithmic method,
tolerances decrease at the slowest rate, resulting in an average reduction of
computation time by 28%. For the linear and exponential methods, tolerances
decrease faster, resulting in speed-ups of only 12% and 0.2%, respectively.
Comparing objective values for the different methods assures that the
reduction of computation time does not compromise the solution quality.

| interpolation | 2 scenarios | 4 scenarios | 6 scenarios
---|---|---|---|---
level bundle | log | 64.6 | 104.1 | 97.3
lin | 166.0 | 81.3 | 111.1
exp | 156.5 | 104.7 | 115.2
none | 136.4 | 109.6 | 111.8
trust-region | log | 69.1 | 71.9 | 67.0
lin | 92.4 | 97.1 | 86.1
exp | 87.7 | 77.0 | 119.2
none | 116.2 | 88.3 | 99.8
Table 1: Impact of interpolation method for inexact cuts on solution time (in
minutes)

Next, Table 2 shows the impact on computation time resulting from the addition
of the valid inequalities to the MP, as described in section 2.3.2. Overall,
the algorithm does not benefit from this refinement increasing the average
solution time by 13%. Therefore, it is reasonable to assume that stabilization
is already sufficient to prevent poor solutions of the MP.

|  | 2 scenarios | 4 scenarios | 6 scenarios
---|---|---|---|---
level bundle | without VI | 64.6 | 104.1 | 97.3
with VI | 116.7 | 86.8 | 102.3
trust-region | without VI | 69.1 | 71.9 | 67.0
with VI | 91.3 | 66.5 | 72.6
Table 2: Impact of valid inequalities on solution time (in minutes)

###  5.3 Benchmark with the deterministic equivalent

This section benchmarks the developed BD algorithm with solving the
deterministic equivalent increasing the number of scenarios up to 16.

It should be noted that the significance of this comparison has its limits
since solving the deterministic equivalent with the barrier algorithm and the
developed BD algorithm utilize computational resources differently. In our
case, the deterministic equivalent runs on a single node with 8 cores and 64
GB memory, the most powerful node available to us, while BD utilizes up to 33
nodes with 4 cores and 16 GB memory each. Considering the benchmark is not
universally valid and depends on the computational resources at hand, we avoid
juxtaposing solution times and rather focus on the scaling behavior.

Fig. 10 shows the run-time to generate and solve the problem for the
deterministic equivalent using between 2 and 16 scenarios. Generation and
solve time both increase exponentially with the number of scenarios including
some outliers, albeit with some outliers.

Figure 10: Performance of the deterministic equivalent

On the other hand, Fig. 11 shows the run-time for the same number of scenarios
using the introduced BD algorithm. Thanks to decomposition and
parallelization, neither the model generation nor the solution time increase
with the number of scenarios. The time for model generation is virtually
constant and only computing the heuristic solution to initialize the
stabilization method adds some fluctuations. Running BD and solving the model
is subject to greater fluctuations, but BD generally outperforms the
deterministic equivalent, if the number of scenarios exceeds four. In
addition, the trust-region method consistently outperforms the level-bundle
method by an average of 33%.

Figure 11: Performance of BD with quadratic trust-region

##  6 Conclusion

This paper applied BD to two-stage stochastic problems for energy planning
under climate uncertainty, a key problem for the design of renewable energy
systems. To improve performance, we modify the standard algorithm and adapt
various refinements to the problem’s characteristics—a simple continuous MP,
and few but large SPs. These refinements include inexact cuts and valid
inequalities, but the focus is on stabilization methods to mitigate the
oscillating behavior of BD. Specifically, we compare proximal and level bundle
methods to a quadratic trust-region approach that is theoretically equivalent
but has not been applied to continuous problems in the literature so far.

In a quantitative case-study stabilization accelerates BD by a factor of four
to ten slightly increasing solve time of the MP, but greatly reducing the
number of iterations. With the computational resources at our disposal, our
refined BD outperforms off-the-shelf solvers for the deterministic equivalent
if planning covers more than four weather years, because run-time does not
increase with the number of scenarios thanks to distributed computing. Among
the three investigated stabilization methods, the trust-region approach
demonstrates the best performance, followed by the level bundle method, which
is on average 33% slower.

By using inexact cuts and setting higher convergence tolerances for the SPs in
the initial iterations, we achieved an average reduction of computation time
by 28%. However, our implementation of valid inequalities did not yield any
improvements.

Although the literature does not offer any obvious further enhancements of BD
for planning renewable energy systems under climate uncertainty, there are
still several promising approaches for further improvement. To present them,
we again follow the categorization of enhancements introduced in Rahmaniani .
(2017):

  * •

Solution procedure: Considering most of the run-time is spent on the SPs,
solving approximations instead of the original SPs can further increase
performance, especially in early iterations. Not solving SPs to optimality in
the current algorithm already exploits that valid cuts do not require exact
solutions. Conceivable approximation methods include surrogate models based on
machine learning, geometric interpolation, or using a reduced temporal
resolution (Neumann Brown, 2021; Mazzi ., 2021; Göke Kendziorski, 2021).
Approximations face two challenges though: First, the high dimensionality of
inputs to the SP that corresponds to the number of complicating capacity
variables. Second, the use of asymptotically inexact oracles also requires a
revision of the stabilization method.

  * •

Solution generation: Although heuristic methods are insufficient for adequate
planning, they can greatly narrow the solution space of the MP to improve
convergence. This approach also offers synergies with the selection of
representative weather years for the planning problem in the first place. For
instance, solving the deterministic problem for specific weather years, as
performed in section 4.2, is not only useful for clustering weather years, but
can also derive lower and upper bounds for capacity variables or system costs.

Furthermore, this paper focuses on basic stabilization methods, but the
literature provides various refinements and, thus, presents additional
opportunities for improvement.

  * •

Cut generation: The disproportional impact storage technologies have on the
run-time that was discussed at the end of section 5.1 suggests the algorithm
would benefit from advanced cut generation. However, the size of the SPs still
poses an obstacle to this approach. Therefore, refined approaches that only
solve a modified but not the original SP are most promising (Papadakos, 2008;
Sherali Lunday, 2011).

Apart from further enhancements, various practical applications benefit from
the introduced BD algorithm. First and foremost, it enables more robust
planning of renewable macro-energy systems by including a large number of
weather years, which was the original motivation for this work. How many
weather years to include and how to select these years for adequate system
planning is an open question for future research. In addition, future research
should extend the optimization under uncertainty to the operational stage in
the SPs. Currently, uncertainty is limited to expansion in the MP, and
operation in each scenario assumes perfect foresight, but stochastic dual
dynamic programming could capture uncertainty at the operational stage as well
(Papavasiliou ., 2018).

Finally, our trust-region method benefits energy planning models that are
deterministic and limited to a single weather year, but capture learning
effects using piecewise-linear cost curves (Zeyen ., 2022). In this case,
linearization adds mixed-integer variables to the expansion problem and as a
result, even standard BD outperforms off-the-shelf solvers for a deterministic
problem (Felling ., 2022). Generally, mixed-integer variables in the MP will
increase the advantage over off-the-shelf solvers and could even speed up BD
when reducing the solution space; for instance, when modeling how the
expansion of large-scale technologies, like nuclear plants, is rather a
discrete than a continuous decision.

##  7 Acknowledgments

The research leading to these results has received funding from the German
Federal Ministry for Economic Affairs and Energy via the project ”MODEZEEN”
(grant number FKZ 03EI1019D). A special thanks goes to all Julia developers.
In addition, we would like to express our gratitude to the reviewers for their
insightful feedback on earlier drafts of this paper.

## Supplementary material

The data and execution scripts for the case-study are available here:
https://github.com/leonardgoeke/EuSysMod/releases/tag/bendersDecompositionPaper.

## References

  * Backe . (2022) Backe2022Backe, S., Skar, C., del Granado, PC., Turgut, O. Tomasgard, A.  2022\.  EMPIRE: An open-source model based on multi-horizon programming for energy transition analyses Empire: An open-source model based on multi-horizon programming for energy transition analyses.  SoftwareX17100877.  10.1016/j.softx.2021.100877
  * Benders (1962) Benders1962Benders, J.  1962\.  Partitioning procedures for solving mixed-variables programming problems. Partitioning procedures for solving mixed-variables programming problems.  Numerische Mathematik4238-252.
  * Bertsimas Mundru (2022) Mundru2022Bertsimas, D. Mundru, N.  2022\.  Optimization-Based Scenario Reduction for Data-Driven Two-Stage Stochastic Optimization Optimization-based scenario reduction for data-driven two-stage stochastic optimization.  Operations Research.  10.1287/opre.2022.2265
  * H. Bloomfield . (2020) bloomfield2020Bloomfield, H., Brayshaw, D. Charlton-Perez, A.  2020\.  MERRA2 derived time series of European country-aggregate electricity demand, wind power generation and solar power generation. Merra2 derived time series of european country-aggregate electricity demand, wind power generation and solar power generation.  10.17864/1947.239
  * HC. Bloomfield . (2016) Bloomfield2016Bloomfield, HC., Brayshaw, DJ., Shaffrey, LC., Coker, PJ. Thornton, HE.  2016\.  Quantifying the increasing sensitivity of power systems to climate variability Quantifying the increasing sensitivity of power systems to climate variability.  Environmental Research Letters1112124025.  10.1088/1748-9326/11/12/124025
  * Bonnans . (2006) Bonnans2006Bonnans, JF., Gilbert, JC., Lemaréchal, C. Sagastizábal, CA.  20065\.  Numerical Optimization: Theoretical and Practical Aspects Numerical optimization: Theoretical and practical aspects.  (2 , 13-161).  Springer-Verlag.
  * Brandenberg Stursberg (2021) Brandenberg2021Brandenberg, R. Stursberg, P.  2021\.  Refined cut selection for benders decomposition: applied to network capacity expansion problems Refined cut selection for benders decomposition: applied to network capacity expansion problems.  Mathematical Methods of Operations Research943383-412.  10.1007/s00186-021-00756-
  * Collins . (2018) Collins2018Collins, S., Deane, P., Ó Gallachóir, B., Pfenninger, S. Staffell, I.  2018\.  Impacts of Inter-annual Wind and Solar Variations on the European Power System Impacts of inter-annual wind and solar variations on the european power system.  Joule2102076-2090.  10.1016/j.joule.2018.06.020
  * Conejo . (2006) Conejo2006Conejo, A., Castillo, E., Minguez, R. Garcia-Bertrand, R.  2006\.  Decomposition Techniques in Mathematical Programming: Engineering and Science Applications Decomposition techniques in mathematical programming: Engineering and science applications.  Springer Berlin Heidelberg.
  * Cordeau . (2006) Cordeau2006Cordeau, JF., Pasin, F. Solomon, MM.  2006\.  An integrated model for logistics network design An integrated model for logistics network design.  Annals of Operations Research14459–82.  10.1007/s10479-006-0001-3
  * Creutzig . (2017) Creutzig2017Creutzig, F., Agoston, P., Goldschmidt, JC., Luderer, G., Nemet, GF. Pietzcker, RC.  2017\.  The underestimated potential of solar energy to mitigate climate change The underestimated potential of solar energy to mitigate climate change.  Nature Energy217140.
  * de Oliveira Solodov (2016) Oliveira2016de Oliveira, W. Solodov, M.  2016\.  A doubly stabilized bundle method for nonsmooth convex optimization A doubly stabilized bundle method for nonsmooth convex optimization.  Mathematical Programming156125–159.  10.1007/s10107-015-0873-6
  * de Oliveira Sagastizábal (2014) de_oliveira_level_2014de Oliveira, W. Sagastizábal, C.  2014\.  Level bundle methods for oracles with on-demand accuracy Level bundle methods for oracles with on-demand accuracy.  Optimization Methods and Software2961180-1209.  10.1080/10556788.2013.871282
  * de Oliveira . (2014) de_oliveira_convex_2014de Oliveira, W., Sagastizábal, C. Lemaréchal, C.  2014\.  Convex proximal bundle methods in depth: a unified analysis for inexact oracles Convex proximal bundle methods in depth: a unified analysis for inexact oracles.  Mathematical Programming1481241-277.  10.1007/s10107-014-0809-6
  * Dupačová, J. and Gröwe-Kuska, N. and Römisch, W. (2003) dupa2003Dupačová, J. and Gröwe-Kuska, N. and Römisch, W.  2003\.  Scenario reduction in stochastic programming - An approach using probability metrics Scenario reduction in stochastic programming - an approach using probability metrics.  Mathematical Programming95493–511.  10.1007/s10107-002-0331-0
  * Felling . (2022) Felling2022Felling, T., Levers, O. Fortenbacher, P.  2022\.  Multi-horizon planning of multi-energy systems Multi-horizon planning of multi-energy systems.  Electric Power Systems Research212108509.  10.1016/j.epsr.2022.108509
  * Frangioni (2002) Frangioni2002Frangioni, A.  2002\.  Generalized Bundle Methods Generalized bundle methods.  SIAM Journal on Optimization131117-156.  10.1137/S1052623498342186
  * Frangioni (2020) Frangioni2020Frangioni, A.  2020\.  Standard Bundle Methods: Untrusted Models and Duality Standard bundle methods: Untrusted models and duality.  AM. Bagirov, M. Gaudioso, N. Karmitsa, MM. Mäkelä S. Taheri (), Numerical Nonsmooth Optimization: State of the Art Algorithms Numerical nonsmooth optimization: State of the art algorithms ( 61-116).  Springer International Publishing.  10.1007/978-3-030-34910-3_3
  * Frangioni Gorgone (2014) Frangioni2014Frangioni, A. Gorgone, E.  2014\.  Bundle methods for sum-functions with “easy” components: applications to multicommodity network design Bundle methods for sum-functions with “easy” components: applications to multicommodity network design.  Mathematical Programming145133-161.  10.1007/s10107-013-0642-3
  * Göke (20211) Goeke2020aGöke, L.  20211\.  AnyMOD.jl: A Julia package for creating energy system models Anymod.jl: A julia package for creating energy system models.  SoftwareX16100871.  10.1016/j.softx.2021.100871
  * Göke (20212) Goeke2020bGöke, L.  20212\.  A graph-based formulation for modeling macro-energy systems A graph-based formulation for modeling macro-energy systems.  Applied Energy301117377.  10.1016/j.apenergy.2021.117377
  * Göke Kendziorski (2021) Goeke2021aGöke, L. Kendziorski, M.  2021\.  The adequacy of time-series reduction for renewable energy systems The adequacy of time-series reduction for renewable energy systems.  Energy238121701.  10.1016/j.energy.2021.121701
  * Göke . (2023) Goeke2021Göke, L., Weibezahn, J. von Hirschhausen, C.  2023\.  A collective blueprint, not a crystal ball: How expectations and participation shape long-term energy scenarios A collective blueprint, not a crystal ball: How expectations and participation shape long-term energy scenarios.  Energy Research & Social Science97102957.  0.1016/j.erss.2023.102957
  * Hilbers . (2019) Hilbers2019Hilbers, AP., Brayshaw, DJ. Gandy, A.  2019\.  Importance subsampling: improving power system planning under climate-based uncertainty Importance subsampling: improving power system planning under climate-based uncertainty.  Applied Energy251113114.  10.1016/j.apenergy.2019.04.110
  * Kaut (2021) kaut2021Kaut, M.  2021\.  Scenario generation by selection from historical data Scenario generation by selection from historical data.  Computational Management Science18411-429.  10.1007/s10287-021-00399-4
  * Kelley (1960) Kelley1960Kelley, JE.  1960\.  The Cutting-Plane Method for Solving Convex Programs The cutting-plane method for solving convex programs.  Journal of the Society for Industrial and Applied Mathematics84703-712.  10.1137/0108053
  * Kiwiel (1990) kiwiel_proximity_1990Kiwiel, KC.  1990\.  Proximity control in bundle methods for convex nondifferentiable minimization Proximity control in bundle methods for convex nondifferentiable minimization.  Mathematical Programming46105-122.  10.1007/BF01585731
  * Lemaréchal . (1995) lemarechal_new_1995Lemaréchal, C., Nemirovskii, A. Nesterov, Y.  199507\.  New variants of bundle methods New variants of bundle methods.  Mathematical Programming69111–147.  https://doi.org/10.1007/BF01585555  10.1007/BF01585555
  * Lemaréchal . (1981) Lemarechal1981Lemaréchal, C., Strodiot, JJ. Bihain, A.  1981\.  On a bundle algorithm for nonsmooth optimization On a bundle algorithm for nonsmooth optimization.  Nonlinear Programming4245-282.  10.1016/B978-0-12-468662-5.50015-X
  * Linderoth Wright (2003) Linderoth2003Linderoth, J. Wright, S.  2003\.  Decomposition Algorithms for Stochastic Programming on a Computational Grid Decomposition algorithms for stochastic programming on a computational grid.  Computational Optimization and Applications242207–250.  10.1023/A:1021858008222
  * Lohmann Rebennack (2017) Lohmann2016Lohmann, T. Rebennack, S.  2017\.  Tailored Benders Decomposition for a Long-Term Power Expansion Model with Short-Term Demand Response Tailored benders decomposition for a long-term power expansion model with short-term demand response.  Management Science6362027-2048.  10.1287/mnsc.2015.2420
  * Lombardi . (2022) Lombardi2022Lombardi, F., Rocco, MV., Belussi, L., Danza, L., Magni, C. Colombo, E.  2022\.  Weather-induced variability of country-scale space heating demand under different refurbishment scenarios for residential buildings Weather-induced variability of country-scale space heating demand under different refurbishment scenarios for residential buildings.  Energy239122152.  10.1016/j.energy.2021.122152
  * Magnanti Wong (1981) Magnanti1981Magnanti, TL. Wong, RT.  1981\.  Accelerating Benders Decomposition: Algorithmic Enhancement and Model Selection Criteria Accelerating benders decomposition: Algorithmic enhancement and model selection criteria.  Operations Research293464-484.  10.1287/opre.29.3.464
  * Marsten . (1975) Marsten1975Marsten, RE., Hogan, WW. Blankenship, JW.  1975\.  The Boxstep Method for Large-Scale Optimization The boxstep method for large-scale optimization.  Operations Research233389-601.  10.1287/opre.23.3.389
  * Mazzi . (2021) Mazzi2021Mazzi, N., Grothey, A., McKinnon, K. Sugishita, N.  2021\.  Benders decomposition with adaptive oracles for large scale optimization Benders decomposition with adaptive oracles for large scale optimization.  Mathematical Programming Computation131683–703.  10.1007/s12532-020-00197-0
  * Neumann Brown (2021) Neumann2021Neumann, F. Brown, T.  2021\.  Broad Ranges of Investment Configurations for Renewable Power Systems, Robust to Cost Uncertainty and Near-Optimality. Broad ranges of investment configurations for renewable power systems, robust to cost uncertainty and near-optimality.  arXiv.  10.48550/ARXIV.2111.14443
  * Ohlendorf Schill (2020) Ohlendorf2020Ohlendorf, N. Schill, WP.  2020\.  Frequency and duration of low-wind-power events in Germany Frequency and duration of low-wind-power events in germany.  Environmental Research Letters158084045.  10.1088/1748-9326/ab91e9
  * Papadakos (2008) Papadakos2008Papadakos, N.  2008\.  Practical enhancements to the Magnanti–Wong method Practical enhancements to the magnanti–wong method.  Operations Research Letters364444-449.  10.1016/j.orl.2008.01.005
  * Papavasiliou . (2018) Papavasiliou2018Papavasiliou, A., Mou, Y., Cambier, L. Scieur, D.  2018\.  Application of Stochastic Dual Dynamic Programming to the Real-Time Dispatch of Storage Under Renewable Supply Uncertainty Application of stochastic dual dynamic programming to the real-time dispatch of storage under renewable supply uncertainty.  IEEE Transactions on Sustainable Energy92547-558.  10.1109/TSTE.2017.2748463
  * Pfenninger (2017) Pfenninger2017Pfenninger, S.  2017\.  Dealing with multiple decades of hourly wind and PV time series in energy models: A comparison of methods to reduce time resolution and the planning implications of inter-annual variability Dealing with multiple decades of hourly wind and pv time series in energy models: A comparison of methods to reduce time resolution and the planning implications of inter-annual variability.  Applied Energy1971-13.  10.1016/j.apenergy.2017.03.051
  * Pfenninger Staffell (2016) Pfenninger2016Pfenninger, S. Staffell, I.  2016\.  Long-term patterns of European PV output using 30 years of validated hourly reanalysis and satellite data Long-term patterns of european pv output using 30 years of validated hourly reanalysis and satellite data.  Energy1141251-1265.  10.1016/j.energy.2016.08.060
  * Rahmaniani . (2017) Rahmaniani2017Rahmaniani, R., Crainic, TG., Gendreau, M. Rei, W.  2017\.  The Benders decomposition algorithm: A literature review The benders decomposition algorithm: A literature review.  European Journal of Operational Research2593801-817.  10.1016/j.ejor.2016.12.005
  * Rey Sagastizábal (2002) Rey2002Rey, PA. Sagastizábal, C.  2002\.  Dynamical adjustment of the prox-parameter in bundle methods Dynamical adjustment of the prox-parameter in bundle methods.  Optimization512423-447.  10.1080/02331930290019495
  * Ritchie Roser (2020) owidenergyaRitchie, H. Roser, M.  2020\.  Energy. Energy.  https://ourworldindata.org/energy
  * Ruhnau Qvist (2022) Ruhnau2022Ruhnau, O. Qvist, S.  2022\.  Storage requirements in a 100% renewable electricity system: extreme events and inter-annual variability Storage requirements in a 100% renewable electricity system: extreme events and inter-annual variability.  Environmental Research Letters174044018.  10.1088/1748-9326/ac4dc8
  * Rujeerapaiboon . (2022) rujeerapaiboon2022Rujeerapaiboon, N., Schindler, K., Kuhn, D. Wiesemann, W.  2022\.  Scenario reduction revisited: fundamental limits and guarantees Scenario reduction revisited: fundamental limits and guarantees.  Mathematical Programming191207–242.  10.1007/s10107-018-1269-1
  * Ruszczyński (1986) Ruszczynski1986Ruszczyński, A.  1986\.  A regularized decomposition method for minimizing a sum of polyhedral functions A regularized decomposition method for minimizing a sum of polyhedral functions.  Mathematical Programming35309–333.  10.1007/BF01580883
  * Ruszczyński (2003) Ruszczynski2003Ruszczyński, A.  2003\.  Decomposition Methods Decomposition methods.  Stochastic Programming Stochastic programming ( 10, 141-211).  10.1016/S0927-0507(03)10003-5
  * Santoso . (2005) Santoso2005Santoso, T., Ahmed, S., Goetschalckx, M. Shapiro, A.  2005\.  A stochastic programming approach for supply chain network design under uncertainty A stochastic programming approach for supply chain network design under uncertainty.  European Journal of Operational Research167196-115.  10.1016/j.ejor.2004.01.046
  * Schramm Zowe (1992) Schramm1992Schramm, H. Zowe, J.  1992\.  A Version of the Bundle Idea for Minimizing a Nonsmooth Function: Conceptual Idea, Convergence Analysis, Numerical Results A version of the bundle idea for minimizing a nonsmooth function: Conceptual idea, convergence analysis, numerical results.  SIAM Journal on Optimization2121-152.  10.1137/0802008
  * Sepulveda (2020) Sepulveda2020Sepulveda, N.  2020\.  Decarbonization of Power Systems, Multi-Stage Decision-Making with Policy and Technology Uncertainty Decarbonization of power systems, multi-stage decision-making with policy and technology uncertainty Massachusetts Institute of Technology. Department of Nuclear Science and Engineering.  https://hdl.handle.net/1721.1/127311
  * Shapiro . (2009) Shapiro2009Shapiro, A., Dentcheva, D. Ruszczyński, A.  2009\.  Lectures on Stochastic Programming: Modeling and Theory Lectures on Stochastic Programming: Modeling and Theory.  Society for Industrial and Applied Mathematics.  10.1137/1.9780898718751
  * Sherali Lunday (2011) Sherali2013Sherali, H. Lunday, B.  2011\.  On generating maximal nondominated Benders cuts On generating maximal nondominated benders cuts.  Annals of Operations Research21057–72.  10.1007/s10479-011-0883-6
  * Skar . (2014) Skar2014Skar, C., Doorman, G. Tomasgard, A.  2014\.  Large-scale power system planning using enhanced Benders decomposition Large-scale power system planning using enhanced benders decomposition.  2014 Power Systems Computation Conference 2014 power systems computation conference ( 1-7).  10.1109/PSCC.2014.7038297
  * Teichgraeber . (2021) Teichgraeber2021Teichgraeber, H., Küpper, LE. Brandt, A.  2021\.  Designing reliable future energy systems by iteratively including extreme periods in time-series aggregation Designing reliable future energy systems by iteratively including extreme periods in time-series aggregation.  Applied Energy304117696.  10.1016/j.apenergy.2021.117696
  * van Ackooij de Oliveira, (2015) Ackooij2014van Ackooij, w. de Oliveira,, W.  2015\.  Level bundle methods for constrained convex optimization with various oracles Level bundle methods for constrained convex optimization with various oracles.  Computational Optimization and Applications volume57555–597.  10.1007/s10589-013-9610-3
  * Van Slyke Wets (1969) Slyke1969Van Slyke, RM. Wets, R.  1969\.  L-Shaped Linear Programs with Applications to Optimal Control and Stochastic Programming L-shaped linear programs with applications to optimal control and stochastic programming.  SIAM Journal on Applied Mathematics174638-663.  10.1137/0117061
  * van Ackooij de Oliveira (2014) van_ackooij_level_2014van Ackooij, W. de Oliveira, W.  2014\.  Level bundle methods for constrained convex optimization with various oracles Level bundle methods for constrained convex optimization with various oracles.  Computational Optimization and Applications573555-597.  10.1007/s10589-013-9610-3
  * van Ackooij . (2016) van_ackooij_inexact_2016van Ackooij, W., Frangioni, A. de Oliveira, W.  2016\.  Inexact stabilized Benders’ decomposition approaches with application to chance-constrained problems with finite support Inexact stabilized Benders’ decomposition approaches with application to chance-constrained problems with finite support.  Computational Optimization and Applications653637-669.  10.1007/s10589-016-9851-z
  * Wang . (2013) Wang2013Wang, J., Wang, J., Liu, C. Ruiz, JP.  2013\.  Stochastic unit commitment with sub-hourly dispatch constraints Stochastic unit commitment with sub-hourly dispatch constraints.  Applied Energy105418-422.  10.1016/j.apenergy.2013.01.008
  * Zakeri . (2000) Zakeri2000Zakeri, G., Philpott, AB. Ryan, DM.  2000\.  Inexact Cuts in Benders Decomposition Inexact cuts in benders decomposition.  SIAM Journal on Optimization103643-657.  10.1137/S1052623497318700
  * Zeyen . (2022) Zeyen2022Zeyen, E., Victoria, M. Brown, T.  2022\.  Endogenous learning for green hydrogen in a sector-coupled energy model for Europe. Endogenous learning for green hydrogen in a sector-coupled energy model for europe.  arXiv.  10.48550/ARXIV.2205.11901
  * Zverovich . (2012) Zverovich2012Zverovich, V., Fábián, CI., Eldon, EFD. Mitra, G.  2012\.  A computational study of a solver system for processing two-stage stochastic LPs with enhanced Benders decomposition A computational study of a solver system for processing two-stage stochastic lps with enhanced benders decomposition.  Mathematical Programming Computation volume4211–238.  10.1007/s12532-012-0038-z