Stabilized Benders decomposition for energy planning under climate uncertainty

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 stabilization, specifically comparing established bundle methods to a quadratic trust-region approach for continuous problems. Anextensive computational comparison shows that all stabilization methods can significantly reduce computation time. However, the quadratic trust-region and the non-quadratic box-step method are the most robust and straightforward to implement. When parallelized, the introduced algorithm outperforms the vanilla version of Benders decomposition by a factor of 100. In contrast to off-the-shelf solvers, computation time remains constant when the number of scenarios increases. In conclusion, the algorithm enables robust planning of renewable energy systems with a large number of climatic years. Beyond climate uncertainty, it can make an extensive range of other analyses in energy planning computationally tractable, for instance, endogenous learning and modeling to generate alternatives.


Introduction
The global energy system must undergo a major transformation and replace 86% of primary energy from fossil fuels as of 2019 with renewable resources to achieve the objectives of the Paris Climate Agreement and limit global warming, (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 et al., 2017).Estimates for the global potential of PV alone range from 1,585 to 50,580 EJ-at least triple the primary consumption in 2019.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 sectors.

System planning under climate uncertainty
Techno-economic planning models are critical tools for analyzing the transformation towards renewable energy systems.Models decide on the expansion and operation of technologies based on representing the flow and conversion of energy in future systems.Commonly, models are linear optimization problems that minimize the costs to satisfy an exogenous demand considering different boundary conditions, for instance, emission limits.Typical applications include developing long-term scenarios for the energy system, assessing technologies in a system context, and analyzing energy policy (Göke, Weibezahn, & von Hirschhausen, 2023).
The fluctuating nature of PV and wind generation is challenging for system planning.In conventional energy systems dominated by dispatchable thermal plants, reliable planning only requires a small number of representative time-steps.However, renewable systems require much higher temporal resolution to capture fluctuations of PV and wind.Additionally, planning must consider energy storage as a key option to balance fluctuations (Levin et al., 2023).Both characteristics significantly increase problem size and complexity, quickly rendering the underlying linear optimization problem computationally intractable (Göke & Kendziorski, 2021).As a result, the literature proposes different techniques to reduce model size while representing renewables accurately, for instance, iteratively adjusting the representative time-periods or limiting high resolution to selected parts of the system (Göke, 2021b;Teichgraeber, Küpper, & Brandt, 2021).
Nevertheless, all methods for tractable planning of renewable energy systems can only consider one representative year of climate conditions, but research shows that renewable supply and energy demand vary substantially across years (Pfenninger, 2017;Pfenninger & Staffell, 2016).Based on a sample of 40 historical years, H. C. Bloomfield, Brayshaw, Shaffrey, Coker, and Thornton (2016) investigate effects on the British power system and conclude robust planning should consider at least ten of these years.Ruhnau and Qvist (2022) study the effect of inter-annual variability on storage requirements for a stylized fully renewable German power system.Compared to a representative year, storage requirements double when considering 35 years of climate data.Ohlendorf and Schill (2020) specifically analyze the frequency of low-wind power events in Germany threatening system adequacy, again concluding that planning models should consider more climatic years.
Overall, climate uncertainty and risk-aware planning are essential for the security of supply in renewable energy systems.Existing heuristic approaches to consider climate uncertainty are limited: Repeatedly solving the same planning problem for different climate conditions can indicate the extent of variability but not provide a generally reliable solution (Grochowicz, van Greevenbroek, Benth, & Zeyringer, 2023).Solving a single problem with representative time-steps from a multi-year sample of climate conditions cannot capture the entire range of climate conditions (Hilbers, Brayshaw, & Gandy, 2019).Accordingly, recent publications call for new planning methods that deploy stochastic optimization to consider the uncertainty of climate conditions, especially as climate change further increases uncertainty (Craig et al., 2022;Doss-Gollin, Amonkar, Schmeltzer, & Cohan, 2023;Plaga & Bertsch, 2023).
Considering climate uncertainty in system planning creates a two-stage stochastic problem: The first stage decides on capacity expansion with uncertainty regarding climate conditions.The second stage reveals the climate conditions and decides on the operation of capacities, e.g., the production of a thermal power plant.Stochastic programming uses independent scenarios with distinct probabilities to represent different climate conditions in the second stage.
So far, two-stage energy planning is limited to a few scenarios representing a short representative period of a week at maximum (Backe, Skar, del Granado, Turgut, & Tomasgard, 2022).This approach limits problem size and keeps the model tractable, but it is unsuitable for planning renewable systems.
Here, scenarios must be long and cover at least one year of climate conditions to model seasonal storage.In addition, scenarios must cover various climate conditions to capture fluctuations and enable reliable planning.The resulting two-stage problem with multiple one-year scenarios is a simple but large linear problem that off-the-shelf solvers can hardly solve.

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 and 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.In energy planning under climate uncertainty, BD decomposes the problem into a master-problem (MP) addressing technology expansion and several mutually independent sub-problems (SPs) for each scenario.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, Castillo, Minguez, & Garcia-Bertrand, 2006).
In existing applications of BD in energy planning models, each SP corresponds to a short time-period, in some models as short as an hour (Brandenberg & Stursberg, 2021;Jacobson, Pecci, Sepulveda, Xu, & Jenkins, 2023;Lohmann & Rebennack, 2017;Skar, Doorman, & Tomasgard, 2014).As already discussed, such decomposition is not applicable when modeling renewable systems dependent on energy storage.In this case, each SP must cover at least one year to account for seasonal storage.Combined with high temporal detail, this results in a simple MP and several large SPs.Furthermore, MP and SPs can contain discrete and continuous variables.However, 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, extensive research proposes different enhancements to improve the original BD.In a review, Rahmaniani, Crainic, Gendreau, and Rei (2017) group them into the following four most important groups, with three of them being relevant to this paper: • 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.In our case, the more relevant SPs can deploy distributed parallel computing.Furthermore, not solving SPs to optimality can produce inexact but valid cuts (de Oliviera & Solodov, 2020;Zakeri, Philpott, & Ryan, 2000).Alternatively, only selected SPs are solved in each iteration to reduce the computational load (van Ackooij & Frangioni, 2018).
• 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.Multi-cut reformulations generally outperform aggregated cuts for energy planning problems since the MP is easy (Jacobson et al., 2023).
-Problem-specific heuristics that provide an initial feasible solution for BD to 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 in the optimum, are added to the MP, narrowing down its solution space (Cordeau, Pasin, & Solomon, 2006).For instance, Wang, Wang, Liu, and Ruiz (2013) add the second-stage constraints of one scenario to the MP, improving convergence.
-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 but have their origins in general convex optimization.For instance, Linderoth and Wright (2003) apply the  ∞ -norm to limit the maximum difference across all variables between a new and the current best solution of the MP; a method first introduced in a general optimization context in Marsten, Hogan, and Blankenship (1975).Level and proximal bundle methods originate from nonsmooth convex optimization and add a quadratic term to the objective function for stabilization (Lemaréchal, Strodiot, & Bihain, 1981;Ruszczyński, 1986).
• Cut generation: Magnanti and 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 that improve 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 established enhancements aim at problems with a difficult MP and small SPs, typically combinatorial problems with discrete complicating variables.Rahmaniani et al. (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 size of the individual SPs that makes our problem computationally challenging.

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. 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. 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 structure for the remainder of the paper is as follows: The next section 2 generally introduces the investigated planning problem and a standard implementation of the Benders algorithm.The following section 3 focuses on different stabilization methods and corresponding parameter strategies.Afterwards, section 4 discusses additional refinements.Section 5 introduces a specific planning problem for an extensive computational comparison of algorithm setups in section 6.The final section concludes and suggests further developments.

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 climatic years.The following subsections first introduce the closed formulation of the problem, followed by the basic BD implementation used to solve it.The notation uses lowercase letters for parameters, uppercase letters for variables, and Greek letters for anything specific to BD.

Problem formulation
Eqs. 1a to 1j provide the closed formulation of the two-stage stochastic energy planning problem.For a set of technologies  ∈ , the problem decides on expanding capacities  , over a set of chronological but not necessarily consecutive years  ∈  .All variables of the problem are continuous and non-negative.
According to 1d, the capacity  , available in each year depends on expansion in the current and previous years provided by    .In Eqs.1e to 1h these capacities constrain the operation of technologies at each time-step  ∈  , in each scenario  ∈ , and within each year.For generation technologies   , the capacity constraint limiting the generation  ,,, includes a scenario-dependent capacity factor  ,, reflecting the available share of capacity.For storage technologies,   , two different capacities constrain three different operational variables: Charged energy   ,,, and discharged energy   ,,, are both constrained by the power capacity   , ; the current storage level   ,,, is constrained by the energy capacity   , .Net supply from generation and storage technologies must match the exogenous demand  ,, in each year, time-step, and scenario, as expressed in the energy balance in Eq. 1i.The balance includes a loss-of-load variable  ,, to avoid an infeasible problem.The storage balance in Eq. 1d computes the storage level at time-step  based on the storage level in the previous time-step  − 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  and the sum of operational costs  , across all years  and scenarios  weighted according to the scenario probabilities   .Expansion costs defined in Eq. 1b depend on the expansion variable and the specific expansion costs   , of each technology.Operational costs defined in Eq. 1c depend on the costs associated with loss-of-load    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 5 and used for benchmarking includes additional elements.First, the model includes multiple regions of expansion and operation.The model can build transmission capacity between regions 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, there are 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 (2021b).

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 stage and second stage of the 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 the operation in the second stage.Increasing the number of years extends both stages, but increasing the number of stochastic scenarios extends the operational problem only.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 with expansion corresponds to the MP, while several independent SPs cover operation in the second stage.Consider the temporal problem structure outlined in Fig. 2 for demonstration.The first stage models two steps of capacity expansion with a 10-year interval.Multi-year intervals are a common modeling strategy for representing transformation pathways, aiming to reduce computational complexity.The second stage includes two operational scenarios for each expansion step.These scenarios represent different patterns of renewable supply and energy demand.In this example, these patterns use historical data for the climatic years 2000 ′ and 2008 ′ .These two years are selected at random here for illustrative purposes.For details on the methods for selecting a sample of climatic years, see section 9.1 in the Appendix.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. Each operational problem can be solved independently for a given set of capacities if no complicating constraint links the variables from different SPs.An example of such a constraint is an emission limit not enforced for each SP separately but jointly across different years and scenarios, like a carbon budget.
Fig. 3 provides the structure of the decomposed problem, analogously to Fig. 1.It shows how the approach splits the operational problem 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.
Alg. 1 presents the standard version of the BD algorithm applied in this paper. is an iteration counter;   and   are the objective function's current lower and upper bound, respectively.We use a multi-cut version 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 for linear planning probelms since the MP is simple (Jacobson et al., 2023).
After initializing all variables, the algorithm starts iterating with a convergence check (Step 1), which tests whether the optimality gap is within a pre-defined tolerance.
Step 2 solves 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  , using the estimator  , .After solving, the algorithm obtains the resulting capacity vector  () and expansion costs  () for the current iteration and sets the lower bound   to the current objective of the MP.In the first iteration,  , is still unconstrained, and the trivial optimum for the MP is zero, meaning no capacity expansion at all.
Algorithm 1: Multi-cut benders algorithm In Step 3.1, the algorithm solves each SP with capacities fixed to the previously computed  () , as described by Eqs.3a to 3i.We will discuss the parallelization mentioned in the algorithm in section 4.3.Afterward, the algorithm gets the actual costs of each SP  ()  , and the dual values  () ,, of the constraint fixing capacities in Eq. 3i.In our case, SPs are always feasible since the loss-of-load  ,, serves as a slack variable for unmet demand.
Step 3.2 computes the Benders cuts Ω ,, defined by Eq. 4 and adds them to the MP.Each cut improves the estimator  , based on the exact result of the SP  ()  , and its gradient  () ,, at the point  () , .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 et al., 2006).In the literature, cuts of this form are termed optimality cuts instead of feasibility cuts that the algorithm adds when the SP is infeasible, which cannot occur in our case.Section 4 discusses the bundle management in Step 4 already included in the algorithm here.
Once all SPs are solved, summing expansion costs  and accurate operational costs  () , in Step 5 gives the total costs at  () .If these undercut the current   ,  () is a new best solution, and we adjust   accordingly.Finally, the iteration counter increases by one, and the next iteration start if the convergence test in Step 1 fails.

Stabilization
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 for two-stage problems.
In the next section, we describe how the BD algorithm introduced in the previous section can incorporate stabilization.Afterwards, we introduce three groups of methods for stabilization and strategies to parameterize them.

Stabilized Benders algorithm
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  , () be the optimal objective value of the SP for investment year  and scenario  as a function of complicating variable choice .This function and its weighted average are polyhedral (and thus convex and non-smooth), which allows constructing a cutting plane model (Shapiro, Dentcheva, & Ruszczyński, 2009, Prop. 2.2 & 2.3).
The cutting plane approximation of  , () at iteration  corresponds to the minimum of the estimator  , restricted by the cuts Ω  ′ ,, , as shown in Eq. 5.Each cut, representing a supporting hyperplane of  , , is defined by Eq. 4. This formulation is equivalent to the piece-wise maximum of all linearizations across the domain of , as the convex optimization literature usually defines the cutting plane model.
The sum of approximations over all SPs φ() () can accordingly be written as Eq. 6.
By sending a candidate solution  () to an oracle, which corresponds to solving the SPs in our case, we obtain ∑ ,    , ( () ) and subgradients  () ,, , which are in the subdifferentials of the respective  , ( () ).Subgradients generalize the concept of gradients to non-differentiable functions. is a subgradient of  () if  () ≥  ()+⟨, −⟩ ∀.The subdifferential of  () is the set of all subgradients at .The subdifferential is a singleton if the function is differentiable at .On this basis, we improve the cutting plane model in Eq. 6 by adding new cuts and then solve the MP, min  () + φ() (), to obtain a new candidate solution  (+1) .
Although the proof of convergence for min φ() → min () is simple (Kelley, 1960), cutting plane algorithms suffer from two computational issues.First, the algorithm is inherently unstable since the capacities  in two consecutive iterates can be extreme points of the feasible space arbitrarily far apart.In energy planning, this problem is particularly pronounced because the feasible region is hardly restricted and completely continuous.Even enforcing the optimum as a starting value may not improve the solution time since the algorithm is initially unable to evaluate the quality of the solution.Second, the size of the model approximating (), in the case of BD the number of cuts, can become prohibitively large as the number of iterations  increases and slow the algorithm (Frangioni, 2020).We emphasize that these issues are mutually reinforcing: Oscillation increases the number of iterations, which in turn slows down the convergence of the approximating model φ() ().
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 to prevent oscillation (Bonnans, Gilbert, Lemaréchal, & Sagastizábal, 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.Removing non-binding cuts after a certain number of iterations addresses the second issue, as further discussed in section 4. Additionally, stabilization reduces the overall number of iterations, thereby mitigating problems associated with excessive growth of the bundle as well.
All stabilzation methods in the literature on (non-smooth) convex optimization build on the cutting plane method and aim to mitigate oscillation by keeping the next iterate  () relatively close to a stability center, a reference point in the (feasible) space of complicating variables.Commonly, this point is part of the sequence of previous iterates { () } =1,…,−1 .In theory, various measures can quantify the distance between the stability center and the next iterate, but the literature typically relies on the Euclidean norm.For the proximal bundle method, discussed below, the literature constructs a "poor man's Hessian" approximates valuable second-order information otherwise absent in the polyhedral setting to justify the Euclidean distance (Frangioni, 2020).
In the problem described in Section 2.1, the capacity vector  is implicitly defined by the selected expansion path (cf.Equation 1d).In practical applications, many elements of the capacity vector are fixed constants.Therefore, defining the stability components below in terms of the lower-dimensional vector of expansion decisions  is computationally advantageous.
An update of the stability center to the current iterate is called a serious step, and an iteration that keeps the current stability center is called a null step.Nevertheless, a null step is still beneficial since it improves the cutting plane model locally.The literature offers various methods to distinguish between null and serious steps.For example, a serious step can be a step that improves the objective by a specified value compared to the current stability center.In an alternative definition, the ratio between the actual improvement determined by () and the improvement estimated by φ() must exceed a certain threshold.In our case, stabilization achieves the best results when the algorithm makes a serious step in each iteration that improves the current best solution, even marginally.
Algorithm 2 represents our general implementation of stabilization with BD.The following sections will extend this general form with specific stabilization methods and corresponding parameter strategies.
Step 0.2 of the algorithm produces an initial feasible solution.In this study, we solve a deterministic instance of the model with the most probable scenario.Section 6 goes more in-depth on testing different initialization strategies.
Step 2 obtains the next iterate by solving the stabilized MP, but the algorithm still solves the original MP in each iteration to obtain a lower bound.The variable  counts the number of serious steps and is an input to the Parameter strategies discussed below.Again, to be concise, we show a parallelized version of the algorithm here but discuss the parallelization in section 4.3.
We did not include a minimum descent target since tests showed that the algorithm does not profit from this method.However, replacing the condition of Step 5 with the equation below can easily implement this feature: where   gives the descent target in iteration .We refer to Frangioni (2020) for an introductory discussion of a commonly used Armijo-type rule.

Bundle methods
There are three relevant bundle methods for the stabililization of cutting plane algorithms and each methods features at least one dynamic parameter that governs the restrictiveness of the stabilization.While all presented methods are theoretically equivalent in the sense that there exists a combination of respective parameter choices that results in the exact next iterate from each bundle method (see Bonnans et al. (2006), Thm.10.7 for a proof), their practical performance hinges critically on the definition of a parameter strategy, also referred to as t-strategy (Frangioni, 2002)).Finding good t-strategies is not trivial, and selecting a poor one may even produce inferior results compared to an unstabilized cutting plane algorithm.For instance, an overstabilized MP, will not benefit from the global nature of the cutting plane model and may take more iterations than an unstabilized version (cf.Frangioni (2020), Fig. 3.2).Strategies suggested in the literature are usually comprised of a set of rules that determines the parameter values for the next iteration dependent on a set of variables describing the algorithm's progress.For example, such a rule could prescribe that the dynamic parameter doubles after l null steps and halves after m serious steps.In this simple strategy l and m are hyperparameters subject to calibration.
In this section, we describe the bundle methods and corresponding parameter strategies.Parameter strategies depend on the method and replace Step 6 in the previously outlined algorithm.Finding good parameter strategies is complex and can be highly problem-specific (Frangioni (2020)).We highlight the hyperparameters used in each strategy before testing specific values in section 6.

Proximal-bundle methods
The first bundle method is called proximal bundle method (PBM) (Lemaréchal, 1977).This variant adds a penalty term to the MP's objective function to prevent moving away from the stability center.Two-stage stochastic problems widely use this method, for example in Oliveira, Sagastizábal, and Scheimberg (2011) Algorithm 2: General parallel stabilized Benders algorithm , ,  (0) ,, , compute Ω 0,, add to stabilized MP and MP;  0,, = 0 Step 2 (Obtain new iterate) Solve stabilized MP and obtain  () ,  () and  ()  do in parallel Step or Zverovich et al. (2012).The stabilized MP is a quadratic program and is given by: However, PBM can go beyond the  2 norm.Pessoa, Sadykov, Uchoa, and Vanderbeck (2018) use a generalized linear PBM approach using a composite penalty term consisting of three terms.This approach does not penalize small deviations determined by the  ∞ -norm, but penalizes deviations outside this confidence interval based on the  1 -norm.As the penalty parameter goes to infinity, this approach turns into the Boxstep method we will introduce in section 3.2.3.To keep the complexity of methods comprehensible, we limited our computational study to the standard PBM method.
Parameter strategies The dynamic parameter   in Eq. 7 in PBM scales the penalty in the objective function with larger values decreasing it; and smaller values increasing it.
Following de Oliveira and Solodov (2016), we use two ruled-based parameter strategies first proposed by Kiwiel (1990).We implement the strategy by inserting the steps of Algorithm 3 into Step 6 of Algorithm 2. Both strategies decrease the weight in response to a serious step and increase it in response to a null step (Step 6.3).However, the methods compute the auxiliary parameter    (Step 6.2) differently.For both methods, the auxiliary parameter mimics the Hessian used in a Newton-type algorithm on a Moreau envelope of the true function .Of course, this function is unknown and we can only approximate the the optimal descent.For details, we refer to Frangioni (2020) and Hiriart-Urruty and Lemaréchal (1993), Chapter XV.
In PBM-1, the auxiliary parameter uses the ratio between the achieved descent and the descent predicted by the model (Kiwiel (1990)).PBM-2 derives directly from a quasi-Newton formula and requires some safeguards to ensure that the numerator of the fraction term remains non-negative (Lemaréchal & Sagastizábal, 1997).Both methods require the following hyperparameters: the starting value  1 , the minimum parameter value   and the scaling factor . ) if method is PBM-1

Level-bundle methods
The level bundle method (LBM) is the opposite of PBM and the trust-region methods discussed below (Frangioni, 2020).Instead of maximizing the predicted descent of the model (  + (  )) − ( (+1) + φ() ( (+1) )) subject to a constraint or a penalty on ‖ (+1) −    ‖, LBM determines a minimum descent target and minimizes the deviation from the stability center needed to reach it (Lemaréchal, Nemirovskii, and Nesterov (1995)).The stabilized MP for LBM is given by: The level parameter   sets the maximum MP objective for the next iterate.It determines level sets of the cutting plane algorithm, i.e., decisions for capacity expansion for which the algorithm estimates function values below   .
The level set is empty, and the problem is infeasible if the objective of the MP is strictly above the prescribed level for all feasible values of .This information is still valuable since   is a lower bound to the overall problem in this case.Therefore, one way to update the lower bound is to set it equal to   if an iteration resulted in an empty level set (van Ackooij & de Oliveira,, 2015).Alternatively, one can update the lower bound by solving the original and stabilized MP in each iteration.One potential advantage of LBM over PBM is that the dynamic parameter   is in the same order of magnitude as the objective functions, which facilitates setting it to appropriate values.
As an extension, de Oliveira and Solodov (2016) propose a doubly stabilized bundle method that combines PBM and LBM.This method automatically chooses between proximal and level iterates.The latter occurs when the level set constraint is binding, while proximal iterations find reasonable solutions within a level set instead.The method has the advantage that the strategy to set   can factor in the dual variable of the level constraint, but due to the method's complexity, we did not include it in our subsequent benchmark.
Parameter strategies We devise two different parameter strategies for LBM.The first one strategy, LBM-1, is a straightforward strategy using a fixed weight, , to compute a convex combination of the current upper and lower bounds of the problem.We make this adjustment in every iteration, regardless of whether a serious or a null step preceded the adjustment.Alg. 4 formalizes this simple mechanism based on Frangioni (2020).As discussed above, we compute the lower bound by solving the unstabilized MP in every iteration.The second strategy, LBM-2, originates from Brännlund, Kiwiel, and Lindberg (1995).Two aspects of this strategy are more refined: First, it only sets the level to a convex combination of the upper and lower bounds if the previous iteration performed a serious step ( > 0).For two, we only increase the level bound in a null step if the level constraint has been binding as indicated by a non-negative dual variable of the constraint   .These refinements ensure that level descent parameter    is kept constant for a short sequence of null steps before it is eventually adjusted (cf.de Oliveira and Solodov (2016) for additional reasoning on the approach).The hyperparameters for LBM-2 are  and   .

Trust-region bundle methods
The final bundle method stabilizes the algorithm using a quadratic trust-region (QTR).The next iterate is constrained to be in a Euclidean hypersphere around the current stability center defined by a certain radius.The quadratically constrained stabilized MP is given by: Algorithm 5: Dynamic adjustment of level parameter   for LBM-2 Note that PBM is the Lagrangian relaxation of the QTR (Hiriart-Urruty and Lemaréchal (1993)), highlighting the theoretical equivalence of the methods.
The existing literature does not use QTRs, but commonly applies a linear alternative using the  ∞ -norm named Boxstep method (Marsten et al., 1975).Although it keeps the stabilized MP linear, Frangioni and Gorgone (2014) show it is inferior to other bundle methods.
Parameter strategies The dynamic parameter for the trust-region methods, QTR and Boxstep, is the radius   constraining deviations from the current stability center measured with the  2 -norm and the  ∞ -norm, respectively.
For QTR, we describe the corresponding steps of the strategy in Alg. 6.Rather than directly adjusting the radius itself, we introduce a scaling factor   , that defines the radius   relative to the  1 -norm of the current stability center.As a result, the radius changes in response either when the dynamic parameter changes or when a serious step changes the stability center.
We illustrate the stabilization based on four consecutive iterations for a stylized example in Fig. 4. It includes three unbounded expansion variables, restricted by a spherical trust-region defined by the constraint of Eq. 9.
The current stability center   is indicated in orange, and the solution at the current iteration  (+1) in yellow or red, depending on whether it improves the current best or not.The minimizer of  () + () is shown in green.In the first iteration, the distance between the origin and the initial stability center determines the sphere's radius.The first iteration in (a) results in an objective above the current best, but the trust-region is binding, as determined in Step 6.2, and remains unchanged.As the cutting plane model improves locally, the second iteration (b) results in a serious step.Re-centering the trust-region for the third iteration leads to a radius adjustment due to the new  1 -norm of the stability center.The third iteration in (c) results in a null step, but the trust-region is no longer binding, and the factor   is reduced by factor .The fourth iteration with a smaller trust-region in (d) does not improve the solution.Since a small radius will result in numerical issues, we impose a minimum on the dynamic parameter   ,   .The hyperparameters include the starting radius scaling factor  1 , the minimum scaling factor   , and the reduction factor .In addition, there is a threshold  to check if the trust-region is non-binding since the original and stabilized MP will rarely have the same objective, even if the trust-region is non-binding due to numerical inaccuracies.
Adding a rule for increasing the radius again is also conceivable because a small starting radius scaling factor  1 can result in slow convergence.However, this implies additional hyperparameters, and in practice, it is much easier to correct a small starting factor if the first iterations show slow progress.
For the Boxstep method, we test a range of different box sizes, which we also denote by .We define the box size relative to the values for the stability center, so the current values impact the absolute boundaries imposed by the trust-region, resulting in a noncubical box (Marsten et al. (1975)).In order to avoid problems with values of zero, we also impose a minimum upper limit.

Additional refinements
In addition to stabilization, there are various methods to improve the poor convergence of the basic BD algorithm.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.
One refined already mentioned in Alg. 1 and 2, is deleting non-binding cuts after a predefined number of iterations.The literature on bundle methods refers to such strategies as -strategies (Frangioni, 2002).In the algorithms,  ,, tracks in which iteration a cut was created or binding the last time is initialized based on the number of the current iteration .
Step 4 loops over all previous cuts, updating  ,, and deleting cuts that were not created or binding within the last  iterations.
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.We scale all equations according to the methodology outlined in (Göke, 2021a) to improve numerical properties.If scaling does not suffice, small values of  ()  , in the SPs are rounded off to zero.The same applies to  ()  ,, 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 sections discuss additional refinements to the stabilized BD in Alg. 2, namely, inexact oracles, valid inequalities, and parallelization.

Inexact oracles
A large body of literature addresses bundle methods and oracle that return inexact information of the SPs ( de Oliveira, Sagastizábal, & Lemaréchal, 2014;de Oliviera & Solodov, 2020;Papadakos, 2008a;van Ackooij & de Oliveira,, 2015).Although it is generally possible that oracles cannot obtain exact information, in our application, the oracles deliberately only return inexact information { , (),  ,, } to reduce the computational time of solving the SPs.
In practice, we implement this feature by deactivating the crossover step of the Barrier algorithm.As a result, the SPs will return a feasible interior point instead of the optimal basic solution.In addition, we test an asymptotically exact oracle strategy by decreasing the Barrier convergence tolerance as the optimality gap closes, eventually reaching the desired optimality level (van Ackooij & de Oliveira,, 2015;Zakeri et al., 2000).This strategy faces a trade-off: On the one hand, increasing the tolerance reduces the computation time of the SPs, but on the other hand, it result in less accurate cuts and a potential increase of iterations.Against this background, Fig. 5 presents three different interpolation methods for this process.In all cases, the algorithm starts at a tolerance of 1 −2 and decreases to 1 −8 , but the reduction follows an exponential, linear, or logarithmic curve.Alg. 2 can easily implement this approach without further changes.
Beyond the implemented approaches, on-demand oracle accuracy could improve the algorithm further (de Oliveira & Sagastizábal, 2014).With this method, interaction with the solver enables an assessment if the next iterate from Step 2 of Alg. 2 will result in a serious step.However, pratical implementations in the context of two-stage stochastic problems focus on cases where the use of multi-cuts significantly impact the MP's computation time.In our problem, this is not the case since the number of scenarios is small compared to the constraints in the MP and even with multi-cuts the computation time of the MP is insignificant Fábián (2013); Wolf, Fábián, Koberstein, and Suhl (2014).

Valid inequalities
Valid inequalities (VI) 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 et al., 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 the examined two-stage problem, the MP does not include an energy balance and will frequently underinvest in generation capacity, particularly in early iterations.To address this, Wang et al. (2013) add one scenario's operational constraints, including the energy balance, 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.10a and 10b to the MP for each scenario.
The first equation 10a enforces a yearly energy balance, setting a lower bound for generation that is binding when hourly generation patterns precisely match demand without relying on storage methods that incur energy losses.The second equation 10b 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 aggregated 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 also aggregate demand and generation by region, assuming lossless and unrestricted transmission as a lower bound.

Parallelization
As shown in Alg. 1 and 2, we parallelize solving the SPs.As a result, we can also solve the original MP in every iteration when solving the SPs to obtain a lower bound for the algorithm at no computational costs since the MP solves much faster than any of the SPs.This computation of the lower bounds is common for LBM methods with simple MPs relative to the SPs (de Oliveira & Solodov, 2016).It also enables us to implement the parameter strategy for the quadratic trust region.
The efficiency of parallelization is the observed speed-up divided by the number of distributed computation nodes, which corresponds to the number of SPs in our case.In the ideal case, parallelization has no overhead, and all SPs require the same time to converge, resulting in an efficiency of one.

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.The public repository linked in the Supplementary material section provides All files to re-produce the case-study.
The applied planning model includes two expansion years, 2030 and 2040, focusing on the power sector.The planned system is not subject to any emission constraints in 2030 but must fully decarbonize by 2040 to cover a diverse range of system setups.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 modeled at an hourly resolution for the entire year.Hydrogen uses a daily resolution; oil and gas a yearly resolution.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 long-term storage is critical since it is key for balancing seasonal fluctuations in decarbonized systems.
Wind and PV capacity 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.

Results
The first two sections compare the different stabilization methods and additional refinements.Subsequently, we benchmark selected configurations to an off-the-shelf solver for the deterministic equivalent.This benchmark solves the case study while varying the number of stochastic scenarios from 2 to 16.The MP and each SP run on independent computing cluster nodes, each with four cores and 16 GB of memory.Core models can differ between runs and nodes.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.We re-ran unexpected performance outliers to ensure the robustness of the results.For solving problems, each node deploys the Barrier implementation of Gurobi 10.1 without crossover and the NumericFocus parameter set to zero.For BD, we delete unused cuts after 20 iterations and choose a convergence tolerance  of 0.2%.The algorithm sometimes fails to converge at smaller tolerances due to rounding.All deployed variations of BD are implemented in the version of the AnyMOD.jlmodeling framework linked in the Supplementary material.

Stabilization
Tab. 1 gives an overview of the tested configurations for the bundle methods and strategies introduced in section 3.For LBM-5, PBM-1, and PBM-2, we copied the configurations proposed in de Oliveira and Solodov (2016).For all other methods, we chose the configurations ourselves.The initial benchmark covers all permutations of values in the table below.Accordingly, we tested six configurations for LBM-1, 15 for LBM-2, and 27 for PBM-1 and PBM-2 each.For box-step the test covers 4 configurations; for QTR 7. Since this totals 80, initially, each configuration is only tested once using a model with four different scenarios.These tests do not use any refinements other than the respective stabilization.Fig. 7 compares the algorithm's run time when deploying the different configurations.Colored symbols indicate the median and average performance for each group of strategies.Without any stabilization, the computation time amounts to 4'600 minutes, so all methods can substantially accelerate the algorithm, but effects greatly vary by method and configuration.As expected, stabilization improves performance by substantially reducing the number of iterations.For instance, the QTR with  1 = 0.01 reduces the number of iterations to 30 from 660 without stabilization.The adverse impact on the MP is negligible since it accounts for less than 0.4% of the total run-time in both cases.Stabilization even decreases the average time per MP since quick convergence reduces the number of cuts added to the MP that continuously increase its size.
Overall, LBM-2 and PBM-2 are less performative than the other methods.For these methods and PBM-1, it is also challenging to identify a relation between the configuration and performance.For instance, in the case of PBM-1 and PBM-2, the first and second best configurations of both Parameter strategies do not share any values.In addition, the spread of computation times is significant.For LBM-1, performance generally improves for smaller values of  that tighten the stabilization, but the overall spread of results is small; results for BOX are similar.QTR also gives similar results until the initial radius  1 reaches a threshold of 0.005.At this point, the radius is too tight, constraining the MP too much and slowing convergence.
In the following, we reduce the configurations to the two best for each method according to Tab. 2.
Fig. 8 compares initialization with the most probable scenarios at a daily resolution on the right, to no initialization on the left, and initialization with the most probable scenarios but an hourly resolution in the middle.We run the tests for the selected subset of configurations, varying the number of scenarios from 2 to 8. Again, the results show a large spread; in some instances, initialization significantly impacts a specific configuration's performance.Such interactions are expected since the starting point is an important determinant for the stabilized iteration.However, predicting the effect initialization will have on a specific configuration is impossible.On average, stabilization with the most probable scenario and a daily resolution performs best.Therefore, we use this configuration in the following.

Additional refinements
In the next step, we analyze the impact other refinements have on the algorithm's performance.First, Fig. 9 shows the impact the inexact-cuts strategies introduced in section 4.1 have on performance.As described, the four compared cases only vary the strategy for the convergence tolerance of the SPs.Since all skip the crossover phase of the Barrier algorithm that moves from an interior point to a basic solution, they all use inexact cuts.Without this feature, computation times are unreasonably slow, preventing an efficient benchmark.
The results show no clear benefit of inexact-cut strategies on performance.For specific scenarios and configurations, inexact cuts do reduce computation time.For instance, with six and eight scenarios, the logarithmic strategies achieve a median reduction of computation time by 33% and never significantly increase computation time.However, no clear trend is observed for the other scenario setups and inexactcut strategies.For some configurations, non-constant strategies can even significantly increase computation time.What is particularly surprising is that in some cases, for example, in all tests with four scenarios, the average time per subproblem increases.In theory, inexact cuts reduce SP time but may increase the number of iterations.Due to the inconclusive results, we presume our benchmark using constant convergence tolerances.
Next, we compare the impact of valid inequalities in the MP as introduced in 4.2 on performance.Again, Fig. 10 compares performance for a different number of scenarios and configurations of stabilization.Overall, valid inequalities increase the spread of solution times without a robust overall reduction.Similar to stabilization, the increase of computation time for the MP is negligible, but we do not observe a consistent reduction of iterations.Even for specific stabilization methods, there is no identifiable trend.For instance, with four scenarios and QTR with a radius  1 = 0.05, valid inequalities more than double computation time, but at six scenarios, reduce time by 16%.Therefore, we will not deploy valid inequalities in the subsequent benchmarks.
Finally, Fig. 11 shows the benefit of parallelization.As discussed in section 4.3, parallelization does not strictly change the algorithm itself but distributes the SPs across different nodes to solve them in parallel.If the overhead is negligible, the refinement will strictly decrease computation time but requires additional computational resources.
Across all tests, the efficiency of parallelization varies between 43% and 68% with an average of 57%.The overhead is negligible, and different solve times of the SPs account for almost the entire efficiency loss.To a certain extent, differences in the underlying years can explain the fluctuations of solve times, but mostly, they are random.As the number of scenarios increases, there is no added risk of outliers delaying the algorithm; efficiency stays constant, and the speed-up increases with the number of nodes.

Benchmark with the deterministic equivalent
In the final section, we perform a concluding evaluation of the refined Benders algorithm and compare it to solving the deterministic equivalent with an off-the-shelf solver.Such comparisons are always at risk of cherry-picking and comparing the best algorithm setup determined in an extensive series of tests against a single test with an off-the-shelf solver.
In the first step, we reduce the variations of stabilization methods and configurations from ten to four, indicated by the blue shade in Tab. 2. This selection builds on the extensive testing with two, four, six, and eight scenarios in the previous sections.Apart from that, all four reference cases use daily initializations, constant convergence tolerances of the SPs, no valid inequalities, and parallelization.Overall, the graph shows a massive speed-up of BD around a factor of 100 thanks to stabilization and parallelization.For two scenarios, the cherry-picked setup takes 29.4 minutes, outperforming the deterministic equivalent with 77.3 minutes, but the more robust reference cases still require between 132.2 (BOX) and 339.0 minutes (LBM-1).With eight scenarios, the reference cases take between 38.2 minutes (PBM-1) and 134.2 minutes (LBM-1) on average, outperforming the deterministic equivalent that requires 119 minutes.
However, the comparison is still biased since we selected the reference cases based on their performance with two to eight scenarios.Therefore, Fig. 13 defined reference case "out-of-sample" with up to 16 scenarios and compares their performance to the corresponding deterministic equivalents in Fig. 14.
Computation time for the deterministic equivalent increases exponentially, although several substantial occur for 14 and 16 scenarios.The graph shows a significant fluctuation of results for the selected proximal and level setups, and methods do not consistently outperform the deterministic equivalent at a certain threshold.For QTR and Boxstep, computation times barely fluctuate, and the selected setups reliably outperform the deterministic equivalent.QTR and Boxstep also showed little fluctuations when comparing different parameter configurations in section 6.1, while the proximal method fluctuated substantially.There is no significant difference in solution quality between the refined BD algorithm and the deterministic equivalent.On average, objective values for the deterministic equivalent are 0.5 smaller.This comparison has limitations since the deterministic equivalent and the parallelized algorithm utilize computational resources differently.In our case, the deterministic equivalent runs on a single node with 24 cores and 16 GB memory each, the most potent node available to us, while BD utilizes up to 33 nodes with 4 cores and 16 GB memory each.Beyond computational resources, performance can vary for other model setups beyond the range tested in this paper.Since the benchmark is not universally valid, the focus should not be on solution times but on scaling behavior.

Conclusion
This paper applied BD to two-stage stochastic problems for energy planning under climate uncertainty, a critical problem for designing 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, valid inequalities, and parallelization, but the focus is on stabilization methods to mitigate the oscillating behavior of BD.We test a broad range of bundle methods and configurations, including a quadratic trust-region, which is theoretically equivalent to other quadratic methods but novel for continuous problems.
In a quantitative case-study, all stabilization methods can significantly reduce computation time compared to the standard algorithm.However, the quadratic trust-region and the non-quadratic box-step method have practical advantages.Firstly, these methods have proven to be robust because their performance is less sensitive to the numeric configuration.As a result, they require less fine-tuning to a specific problem to achieve good results.Second, the practical implementation of these methods is more straightforward.
The developed algorithm can also benefit from inexact cuts and valid inequalities, but there are also adverse effectand no clear trend.On the other hand, parallelization consistently improves performance by a factor corresponding to the number of SPs divided by two, which translates to a parallelization efficiency of 50%.Overall, the refined algorithm accelerates the vanilla BD algorithm by a factor of 100.Thanks to parallelization, its computation time remains constant when the number of scenarios increases.Cherry-picked setups can already outperform off-the-shelves solvers for the deterministic equivalent for two scenarios; more general algorithm setups start to achieve this around six scenarios, but note that performance can vary for other model setups beyond the range tested in this paper.
There are further conceivable improvements to the algorithm described in this paper.To list them, we again follow the categorization of enhancements introduced in Rahmaniani et al. ( 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.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 (Göke & Kendziorski, 2021;Mazzi, Grothey, McKinnon, & Sugishita, 2021;Neumann & Brown, 2021).However, approximations face two challenges: First, the high dimensionality of inputs to the SP corresponds to the number of complicating capacity variables.Second, using asymptotically inexact oracles also requires a revision of the stabilization method.The use of on-demand accuracy oracles appears promising.Introducing a solver component that allows to skip iterations and replace by them with approximations if they do no meet a predefined descent target is a natural extension of the current algorithm (van Ackooij & de Oliveira,, 2015).
• 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 climatic years for the planning problem in the first place.For instance, solving the deterministic problem for specific climatic years, as performed in section 9.1, is not only useful for clustering climatic years, but can also derive lower and upper bounds for capacity variables or system costs.
• Cut generation: The disproportional impact storage technologies have on the run-time that was discussed at the end of section 6.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, 2008b;Sherali & Lunday, 2011).
Apart from further enhancements, the outlined algorithm enables robust planning of renewable macroenergy systems wit a large number of climatic years, which was the original motivation for this work.How many climatic 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, Mou, Cambier, & Scieur, 2018).
Beyond from climate uncertainty, the algorithm has the potential to make a broad range of other analyses in energy planning computationally tractable.In contrast to the deterministic equivalent, the algorithm scales well when adding complexity to the expansion problem since its computation time is negligible in the decomposed algorithm.For instance, the algorithm is suited for planning with endogenous learning, which adds mixed-integer variables to the expansion problem for a piecewise linear approximation of cost curves (Felling, Levers, & Fortenbacher, 2022;Zeyen, Victoria, & Brown, 2023).Similarly, the algorithm makes a stochastic representation of investment costs or discrete investment decisions in highly resolved models viable.Furthermore, it can efficiently compute near-optimal solutions due to its iterative nature and since oracles remain valid when changing the objective function.

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.data for 39 climatic years from 1980 to 2018 published in H. Bloomfield, Brayshaw, and Charlton-Perez (2020), based on re-analysis data from the NASA's MERRA-2 dataset.
Stochastic programming often reduces  scenarios of historical data, in our case 39 climatic years, to a smaller number  <  that captures the empirical distribution but improves computational tractability (Dupačová, J. and Gröwe-Kuska, N. and Römisch, W., 2003;Rujeerapaiboon, Schindler, Kuhn, & Wiesemann, 2022).Scenario reduction methods determine a set of representative scenarios  and a set of corresponding weights {  ∶  ∈ , ∑ ∈   = 1,   ≥ 0}, such that the weighted scenario average of second stage costs approximates the sample average over all  scenarios.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  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 6 below, it is sufficient to ensure the  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 and Mundru (2022), for selecting  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 ℝ + , the space of system costs.To perform the clustering, we calculate a matrix of symmetric distances between scenarios, or climatic years, where the -th row and  * -th column is computed according to Eq. 11: where  , * refers to total system costs when the deterministic problem for scenario  is solved with capacities fixed to results of the deterministic problem for scenario  * .For instance,  1980,2008 corresponds to total system costs, including infeasibility costs, when solving the deterministic model for the climatic year 1980 with capacities originally computed for 2008.The clustering algorithm selects  representative scenarios and assigns each of the  −  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. 15 illustrates the results of the clustering algorithm for  = 6.Each color represents a group and the node of the representative scenario for this group is equipped with its weight   , which reflects the relative group size.The color and thickness of connecting lines between the nodes indicate similarity, the inverse of the distance  , * as defined above.

Figure 1 :
Figure 1: Structure of closed two-stage problem

Figure 2 :Figure 3 :
Figure 2: Graph of years and scenarios for exemplary problem

Figure 5 :
Figure 5: Interpolation methods for convergence tolerance of SPs

Figure 6 :
Figure 6: Graph of model elements

Figure 7 :
Figure 7: Comparison of configurations for stabilization methods

Figure 8 :
Figure 8: Comparison of initialization strategies

Figure 10 :
Figure 10: Impact of valid inequalities

Figure 12 :
Figure 12: Overall benchmark of refined algorithm

Figure 14 :
Figure 14: Performance of deterministic equivalent and variable costs of generation technologies   ,, , for instance, fuel costs.

Table 2
Two best configurations for each stabilization method and Parameter strategies