Applying High-Performance Computing to the European Resource Adequacy Assessment

This work considers the European Resource Adequacy Assessment, which is a pan-European resource adequacy process that is being developed by the European Networks of Transmission System Operators for Electricity (ENTSO-E). A critical part of this process is the so-called Economic Viability Assessment model, which aims at determining future expansion and retirement capacity opportunities for the entire European network. As such, the problem is stochastic. Nevertheless, due to computational constraints, simplified approaches have been followed by ENTSO-E. Our work formulates the problem as a two-stage stochastic problem and proposes two decomposition algorithms for solving the problem which are implemented in a high-performance computing infrastructure. The first is a subgradient-based algorithm, and the second uses a relaxation of the second stage (the economic dispatch) in order to speed up the subgradient calculation thus achieving a considerable reduction in solution time. We compare our schemes against the commonly used Bender's decomposition. We compare the obtained stochastic solution against the deterministic solution proposed by ENTSO-E for their 2021 study and analyze the impact of the stochastic solution on various adequacy indicators.


I. INTRODUCTION
W ITHIN an uncertain world, measuring and analyzing the ability of the electric power system to react to adverse uncertain conditions has become increasingly important.The Clean Energy Package [1] has recognized the importance of this task with Regulation (EU) 941\2019 [2] and Regulation (EU) 943\2019 [3].The latter stipulates the need for a robust European resource adequacy assessment that provides an instrument for detecting and measuring adequacy concerns [4].In particular, resource adequacy concerns identified through the European Resource Adequacy Assessment (ERAA) are to become the basis for justifying the implementation of capacity mechanisms within European Member States.As required in Regulation (EU) 943\2019 [3], Member States wishing to introduce capacity mechanisms must do so on the basis of an adequacy concern that is identified in the ERAA study, complemented possibly with a national resource adequacy study.Consequently, there is an institutional urge to develop a reliable and robust ERAA study at a pan-European level.
The European Network of Transmission System Operators for Electricity (ENTSO-E) is the body mandated by regulation to develop the methodology and conduct the study on an annual basis [5].The ERAA is a pan-European resource adequacy assessment for up to 10 years ahead, which aims at measuring the ability of the power system to react to a set of future uncertain events [4].The ERAA study covers the entire pan-European interconnected system, thus modeling 56 bidding zones in 37 countries.Adequacy concerns are identified with the help of adequacy indicators, with the Loss of Load Expectation (LOLE) being the most common indicator used among the EU Member States to define the respective reliability standards.The Expected Energy not Served (EENS) indicator is also computed in the ERAA study, in order to assess the depth (MW) of the curtailments.Such indicators are naturally linked to the installed capacity mix.For the purpose of determining the installed capacity mix, the ERAA study introduces the so-called Economic Viability Assessment (EVA), which aims at modeling economic parameters that affect the available generation capacity within Europe.Unfortunately, due to the scale of the ERAA, incorporating an EVA that integrates the stochastic nature of future events has thus far been considered as being out of scope.Therefore, a deterministic approach has been followed in ERAA 2021 [6], while in ERAA 2022 a tractable formulation is obtained by considering a reduced stochastic model, this simplification consists of using 3 of 35 uncertainty realizations.ENTSO-E is moving towards a stochastic programming formulation [7], endorsed by ACER, but currently limited to three scenarios.This motivates the development of a framework which allows tackling the EVA in its stochastic version.
Our work aims at bridging this methodological gap by proposing a novel parallel computing algorithm, based on ideas from stochastic programming, and implemented on high-performance computing infrastructure.Our approach allows us to account for the stochastic nature of the EVA study in ERAA.Furthermore, we study the potential impact of ignoring the stochastic nature of the EVA on the capacity mix, and in turn, the consequences that this could have on adequacy indicators.
The EVA aims at determining capacity expansion and capacity retirement opportunities for the entire European network.As such, it relates to two streams of literature: (i) stochastic capacity generation expansion, with the added complexity of considering retirement opportunities; (ii) large-scale optimization, which is tackled with the aid of high-performance computing.

A. Stochastic Capacity Generation Expansion
The stochastic capacity generation expansion literature presents a variety of strategies for addressing generation expansion. 1When the size of the problem is manageable, the problem is solved directly by a commercial solver [8], [9], [10], [11].In order to decrease the computational burden, certain authors have followed scenario selection techniques, and report solving time improvements with corresponding losses in the quality of the solution [10], [12], [13], [14], [15], [16], [17].Note that these approaches typically still rely entirely on a commercial solver for solving the reduced problem.By relying on stochastic programming formulations [18], certain authors use Bender's decomposition in order to tackle the problem [19], [20], [21], [22].Furthermore, modelling tools have been developed in order to solve the problem pragmatically.Examples of open-source software include the EMPIRE model [23] which solves the problem as a large-scale linear problem, while [24] has proposed extensions of EMPIRE where the problem is decomposed by using progressive hedging.Recent approaches [25], applied to the Brazilian power system, have proposed extensions to Benders decomposition in order to speed up convergence.Examples of commercial tools include the PLEXOS software, which allows solving the problem as a large-scale linear problem or using Bender's decomposition.Certain heuristics have also been proposed, including a rolling-horizon scheme [26] and an hourly aggregation of time series [27].
The modelling specifications of the EVA prevent us from solving the problem directly using a commercial solver.These specifications result in a large-scale problem, due to its wide geographical scope, combined with the time step chosen by ENTSO-E, as well as the endogenous representation of uncertainty.We highlight that a small set of uncertainty realizations are already enough to produce a problem of considerable size, thus preventing the use of a scenario reduction technique which allows tackling the problem in its extended form without using a decomposition strategy.On the other hand, a drawback of Bender's decomposition is that its performance diminishes as the number of expansion/retirement candidates increases.The progressive hedging algorithm presented in [24] appears to deteriorate in performance as the size of each scenario subproblem increases, which renders it impractical for our purposes.Some early theoretical work from 1986 [28] presented the possibility of employing subgradient schemes in deterministic capacity expansion problems.Furthermore, past research has proven the effectiveness of subgradient algorithms when dealing with large instances of stochastic unit-commitment [29], [30].Consequently, our work aims at translating this success into the ERAA study.Due to the EVA modelling specifications, the calculation of each subgradient is time-consuming.Therefore, our work, in addition to putting forward a subgradient-based algorithmic framework, proposes a further algorithm that calculates approximations of the subgradients efficiently, thus providing notable computational benefits.

B. High Performance Computing
High-performance computing has become critical in tackling large-scale power system optimization problems.Early work on the topic [31], [32] describes parallel computing schemes for addressing security-constrained optimal power flow and hydro systems, respectively, using Bender's decomposition.Parallel computing schemes for Lagrangian decomposition have been developed for solving optimal power flow problems [33], [34].In recent years, parallel computing has enabled tackling large-scale instances of stochastic unit commitment [35].In particular, the usage of asynchronous parallel computing has resulted in a reduction of computation time from weeks to a few hours for certain stochastic unit commitment instances [30].The hydrothermal scheduling literature has benefited from the usage of high-performance computing as well.While the known Stochastic Dual Dynamic Programming (SDDP) algorithm [36] has been developed for tackling such a problem, as explained in [37] SDDP can encounter challenges, and several parallel computing studies have followed.For instance, in [38] the authors propose a synchronous parallel scheme, while [39], [40] have proposed asynchronous implementations.In [41] the authors have compared these different parallel implementations.Some recent work [42] has shown further parallelization strategies that are capable of outperforming the commercial implementation of PSR-SDDP 2 which is targeted for hydrothermal scheduling.
The parallelization of decomposition techniques such as Bender's decomposition or subgradient methods has been studied in the past [31], [32], [34], [35], [43].Nevertheless, the stochastic capacity expansion literature has been short in exploiting the benefits of parallel computing.In this work, we bridge this 2 PSR is a consulting firm based in Rio de Janeiro that has pioneered the commercialization of the SDDP algorithm for hydro-thermal planning methodological gap by proposing a parallelization strategy for the considered algorithms, which allows us to arrive at solutions for the stochastic formulation of EVA within a few hours of computation.

C. Organization & Contributions
(i) In terms of methodological contribution, we solve the stochastic EVA considering all available uncertainty conditions.(ii) Our algorithmic contribution amounts to proposing a parallel computing subgradient-based method and a novel parallel computing second-stage relaxation scheme, which decreases computational burden significantly.These algorithms are benchmarked against the commonly used Benders decomposition (which is also the algorithm that the PLEXOS commercial software uses).(iii) Our policy analysis contribution amounts to further testing the quality of the solution by comparing it against the ERAA 2021 deterministic approach.
The organization of the work is as follows.The EVA is formulated as a two-stage stochastic capacity expansion problem in Section III.This formulation leads to a large-scale stochastic programming problem.A customized algorithm for tackling such a problem is proposed in Section IV.Section V proposes parallel schemes for the described algorithms.Finally, the results are discussed in Section VI.

II. EUROPEAN RESOURCE ADEQUACY ASSESSMENT
The European resource adequacy assessment proposes a methodology for measuring the ability of the power system to react to uncertain events [44].The overall methodology consists of two main blocks.The first block aims at determining investment and retirement opportunities (which we refer to as the expansion plan), the so-called Economic Viability Assessment (EVA).The second block uses these opportunities in order to measure adequacy indicators.
Fig. 1 provides a visual representation of the overall ERAA methodology.During the first step (the EVA), a variety of uncertain climate conditions are considered (ENTSO-E considers a total of 35 climatic years for ERAA 2021 and ERAA 2022).The objective is to decide on investment and retirement that minimize the expected operational cost of the system.Note that this has to be decided before the realization of uncertainty, thus leading to a stochastic problem.During the second step, the adequacy indicators are calculated.The previously computed expansion plan is used for this purpose, and it is evaluated over an uncertainty set that consists of climate years as well as random outage patterns.
We highlight a critical difference between these steps.Step 2 can be decoupled into several independent problems, one for each climate year and outage pattern.This implies that it is computationally tractable.The situation regarding step 1 is considerably more involved, as it naturally links all climate years into a single problem, meaning that it is not possible to decouple the climatic years into independent problems.In our work, we focus on the first step.In order to tackle it, we propose a parallel computing algorithmic framework.Due to increased computational complexity, ENTSO-E has considered simplified approaches in order to tackle step 1.For ERAA 2021 the simplification has been two-fold.(i) On the one hand, a scenario reduction methodology is applied which selects 7 out of 35 climatic years.(ii) The optimal expansion plan for each one of the selected climatic years is calculated, thus obtaining 7 expansion plans.The average between these expansion plans is selected as the approximate solution of the stochastic problem.For ERAA 2022, a stochastic model is proposed.However, due to increased computational complexity, the model is reduced to 3 out of 35 climatic years, which are selected using a scenario reduction technique.
ENTSO-E has provided access to the ERAA 2021 input data for our team in the context of this work.For ERAA 2021, the EVA is modeled as a two-stage problem: the expansion plan is decided for a single target year.Note that ERAA 2022 adopts a multi-stage approach, namely there are consecutive years and the expansion is decided just before each year.In this work, we focus on ERAA 2021, which is the more natural step, since not even this model can be tackled in its stochastic form by state-of-the-art solvers.
The ERAA 2021 data considers the entire interconnected pan-European network, which is represented by considering 56 bidding zones that correspond to 36 countries.This is depicted in Fig. 2. In terms of generators, the data contains the installed capacity mix per zone.We highlight that the generator data is aggregated per technology, per zone.A detailed overview of the installed capacity mix per zone can be found in [45].The generator data also contains time series of planned maintenance of generators, de-rating factors, and must-run profiles of the generators.The data represents the transmission network as a transportation network, and includes transmission lines between  zones, together with their limits.The data contains several candidate retirement opportunities, per zone, for thermal generators.Two thermal investment candidates, per zone, are considered.The uncertainty arises in the form of the so-called climatic years, each one being a time series (per zone) of demand, PV, wind, and inflow profiles.ENTSO-E has considered 35 climatic years for the EVA of the ERAA 2021 edition, which are depicted in Fig. 3.
ENTSO-E uses the ERAA 2021 data as input for their modeling tool for compiling the EVA expansion problem.As this tool is unable to tackle the problem in its stochastic form, we instead use the exact same input data in order to put together an open-source Julia [46] version, which enables us to develop a decomposition scheme for solving the problem, using the full set of 35 climatic years.We highlight that both models are built with the same modelling specifications.In particular, similarly to the EVA ENTSO-E model, ours is a linear model which takes into account the installed capacity mix per zone for proposing investment and retirement decisions (meaning that it is not a greenfield model).

III. EXPANSION PROBLEM
The stochastic capacity expansion problem is formulated as a two-stage stochastic program [18].The first stage determines investments and retirements of technologies.The second stage solves an economic dispatch over a target year.The introduction of uncertainty takes place during the second stage.Each uncertainty realization corresponds to a so-called climatic year, which consists of a time series of demand, solar production, wind production, and hydro inflows.In order to ease the exposition, the sets, variables, and parameters of the problem are presented in the appendix.Furthermore, parameters are denoted in upper case while optimization variables are in lower case.

A. Expansion Constraints
There is a maximum amount of plausible invested/retired capacity.This is modeled with the following constraints:

B. Generator Constraints
Given ω ∈ Ω, the minimum and maximum power generation capabilities of units are described by the following constraints: Constraints 3, 4 model the power production of new and existing capacity.Constraint 5 models must-run obligations.

C. Transmission Network
The transmission network is modeled as a transportation network.The constraints are then bounds on the transfer capacity.Given ω ∈ Ω, these can be described by the following constraints: Transportation models fail to accurately represent the true physics of power flow in European network models.This has motivated the introduction of a zonal PTDF approximation in European market clearing models.This implies the addition of linear constraints, which do not affect the overall structure of the decomposition schemes that are proposed in this article.

D. Batteries
Batteries are modeled as energy storage resources available per zone.For a given ω ∈ Ω, they are modeled as follows: bd n,t,ω ≤ BD for all n ∈ N , t ∈ T (10) Constraint 7 models the load balance of the battery.Constraints 8, 9, 10 model the maximum capacity, charge, and discharge capabilities of the battery.

E. Hydro Power
Hydropower generation is modeled using four different hydro technologies: run-of-river, reservoir, pumped storage open loop, and pumped storage closed loop.The hydrology is simplified by ENTSO-E by considering the aggregated hydrological capabilities of each zone.The inflows are measured as equivalent energy inflows (MW).For a given ω ∈ Ω, each one of these technologies is modeled as follows.For ease of exposition, a variable associated with a certain technology has an associated subscript: run-of-river R, reservoir S, pumped storage open loop O, and pumped storage closed loop C.
r Run-of-River: There is no storage capability, therefore the net inflows are considered as turbined water.
r Reservoir: There is storage capability, consequently the water can be turbined or stored.
r Pumped storage open loop: There are head and tail reser- voirs.The head turbines water to the tail reservoir, thus producing power.In low-demand periods the tail pumps water back to the head reservoir, thus consuming power.
The system is exposed to natural inflows.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Constraints 15, 16 are the water balance constraints of the head and tail reservoirs respectively.Note that the head reservoir is subject to rainfall uncertainty.Constraint 17 bounds the maximum storage, while constraints 18, 19 bound the turbined and pumped water respectively.
r Pumped storage closed loop: this is modeled in the same way as the open loop case, with the difference that no natural inflows are considered.
for all n ∈ N , t ∈ T (24)

F. Load Balance
The load balance constraint aims at satisfying the demand of each zone using the resources of each zone plus imports from neighbouring zones.For a given ω ∈ Ω, it is formulated as follows: We highlight that reserve requirements are modelled in the ERAA 2021 study as extra load.More accurate models for reserve requirements have been studied in the past [47], their inclusion implies the addition of linear constraints which do not disrupt the overall structure of the algorithms developed in this article.

G. Objective Function
For a given ω ∈ Ω, we define the following quantities: n∈N ,g∈G,t∈T n∈N ,g∈G * ,t∈T n∈N ,l∈L(n),t∈T Equation ( 26) is the first-stage cost.It consists of the investment cost plus the fixed maintenance cost of new capacity x * n,g , minus the fixed maintenance cost of retired capacity x n,g .Equation ( 27) corresponds to the cost of producing p n,g,t,ω units of power, similarly (28) corresponds to the generation cost of new capacity.Equation ( 29) corresponds to the cost of transporting power through the transmission network.Equation (30) corresponds to a penalty for water spillage.Finally, (31) is the cost of involuntary load shedding, which is penalized at VOLL.Putting these elements together leads to the stochastic capacity expansion problem, which is described as follows: 28)+( 29)+( 30)+( 31) (3)-( 25) for all ω ∈ Ω (CEP) Note that we have distinguished between two sets of constraints.The former refers to the so-called first-stage constraints, therefore they do not depend on ω ∈ Ω.The latter set of constraints is the second-stage constraints, they do depend on uncertainty, and there is one such set of constraints for each ω ∈ Ω.Note that this formulation implies that the first-stage variables x n,g , x * n,g are decided before uncertainty realizes, and thus do not depend on ω ∈ Ω.

IV. SOLUTION STRATEGY
The stochastic capacity generation expansion literature has often relied on Bender's decomposition (also known as the L-Shaped method when uncertainty is introduced) [18].This scheme performs poorly as the number of expansion/retirement possibilities increases.In view of this, we propose a subgradient algorithm that is better suited for such situations.This technique allows us to reduce the number of iterations.However, it can still be computationally costly as the calculation of each subgradient is non-trivial.Consequently, we further propose a relaxation of the economic dispatch, which allows us to calculate subgradient approximations efficiently.The approximation is refined throughout iterations, thus ensuring convergence.Such a technique allows us to reduce the computational burden significantly.We begin by briefly describing the L-shaped technique, followed by our algorithmic contributions.
The L-shaped method breaks down the overall problem into smaller subproblems, first by considering separate subproblems for the first and second stages, and secondly by considering a subproblem for each uncertainty realization for the second stage, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.as presented in Fig. 4. Given an uncertainty realization ω and a first-stage decision x, x * , the second-stage subproblem for a given ω is as follows: 28) + ( 29) + ( 30) + ( 31) where λ * n,g,t,ω , λ n,g,t,ω are dual multipliers of constraints ( 3), (4) respectively.These subproblems allow us to re-write the CEP problem as The function E ω [V(x, x * , ω)] is piece-wise linear convex in x, x * [18], and can therefore be under-approximated by supporting hyperplanes, commonly referred as to cuts.These cuts are computed by calculating the subgradients of the second-stage functions V(x, x * , ω) [18].Consequently, given a collection of supporting hyperplanes {c i (x, x * )} N i=1 , the previous problem can be approximated as follows: The L-Shaped algorithm advances by finding, at each iteration, a new supporting hyperplane for problem M. Each iteration begins by solving problem M which leads to a trial action x i , x * i .The second-stage subproblems V(x, x * , ω) are then solved around x i , x * i .The dual multipliers λ n,g,t,ω , λ * n,g,t,ω are a subgradient of V(x, x * , ω) at x i , x * i , and can therefore be used for computing a supporting hyperplane [18], which is added to problem M. As we move through iterations, the method starts building an accurate representation of E ω [V(x, x * , ω)] around the optimal region, eventually finding the optimal value.In fact, the method converges after finitely many iterations [48].

A. Subgradient Algorithm
A drawback of the L-Shaped scheme is that, as the dimensionality of x, x * increases, additional supporting hyperplanes are required in order to describe E ω [V(x, x * , ω)].As a consequence, finding the optimal region may require many extra costly iterations.Instead, in this work we use a subgradient algorithm.Given an initial expansion candidate x, x * , we specifically calculate a subgradient of the objective function of the CEP-R problem around x, x * and update the candidate expansion plan along the direction of such a subgradient.This method has the following advantages: (i) it does not require a hyperplane description of E ω [V(x, x * , ω)] in order to advance to the next candidate x, x * ; (ii) it can be initialized around a trial x, x * which is known in advance to be somewhat close to the optimal value, thus ensuring that if the starting value is close to the optimal solution the iterates remain in the optimal region and will need few iterations.
We start by decomposing problem CEP by rewriting it as CEP-R.Note that a subgradient of the objective function along the x * n,g coordinate is: . And a subgradient along the x n,g coordinate is: Due to the reformulation CEP-R, the calculation of these slopes can be decomposed into calculating several subproblems, concretely by solving V(x, x * , ω) for each ω ∈ Ω.This decomposition allows us to apply the following scheme.We commence by providing an initial candidate action x 1 , x * 1 .During each iteration, the subproblems V( x i , x * i , ω) are solved for all ω ∈ Ω.Using the dual multipliers of these subproblems, the subgradients ρ n,g , ρ * n,g are calculated.Finally, the trial action is updated through a projected subgradient step: The term α i is a stepsize that is crucial to the performance of the algorithm [49].We have selected the Polyak stepsize [49], which ensures convergence.The Polyak stepsize can be described as follows: Here, W * is the optimal value of the CEP problem, while W i is the current iterate objective value.As the optimal value W * is not known in advance, an approximate value can be used.No lower bound is obtained, thus its stabilization is used as a stopping criterion.

B. Second-Stage Relaxation Algorithm
Each iteration of the subgradient scheme can be costly because it carries the computational burden of solving each second-stage subproblem.For this reason, we propose a scheme that relaxes the second stage, namely the economic dispatch.Using such a relaxation, each iterate of a trial x i , x * i becomes efficient, thus allowing us to increase the search speed.The relaxation can then be refined in order to tighten the search and ensure convergence.
Let us begin by describing the second-stage relaxation.To achieve this, we resort to dynamic programming [50], and partition the second-stage horizon 1, . . ., T into K consecutive blocks: {1, t 1 }, {t 1 + 1, t 2 }, . . ., {t K−1 + 1, T }.This leads to the representation shown in Fig. 5, where the second stage has been partitioned into several blocks.The subproblems at each block are given by the dynamic programming equations.The equation at block k is given by: 27) + ( 28) + ( 29) + ( 30) + ( 31 Here, the notation indicates that the objective function and the constraints are restricted to t = t k−1 + 1, . . ., t k .Note that, at the initial time boundary, i.e. at t = t k−1 + 1, constraint (7) requires the battery state of charge at t = t k−1 .Consequently, the notation indicates that V k depends on bv t k−1 ,ω , the battery state of charge at stage t k−1 .Similarly, ( 12), ( 15), ( 16), (20), (21) require the reservoir level at stage t k−1 , and this information is summarized in the vector v t k−1 ,ω .At the final time boundary at t = t k , the objective function includes V k+1 , which captures the future costs of the system.Note that the subproblem of the first block satisfies Each function V k+1 is piecewise linear convex in x, x * , bv t k ,ω , v t k ,ω [18], and so we can approximate it using supporting hyperplanes.Thus, given a collection of supporting hyperplanes {c i (x, x * , bv t k ,ω , v t k ,ω )} N i=1 an approximation of the subproblem at the k-th block is given by: 27) + ( 28) + ( 29) + ( 30) + ( 31 In particular, the approximation of the first-block subproblem V 1 is an approximation of the second-stage subproblem, namely the economic dispatch.Consequently, we can approximate problem CEP-R as follows: Calculating each V 1 is straightforward, therefore problem CEP-A can be solved efficiently.Note, however, that, due to the approximation, the subgradients we may obtain are not necessarily tight.Consequently, the approximation is tightened throughout iterations, thus ensuring convergence.An initial approximation of V 1 (x, x * , ω) is calculated at the beginning of the algorithm.The objective of doing so is to, similarly to the subgradient scheme, have an initial candidate expansion plan to start the search.These ideas are the basis of our algorithmic scheme, which is depicted in pseudo-code in Algorithm 1.
The second-stage relaxation algorithm is illustrated graphically in Fig. 6.The figure presents step 2) of the algorithm.For ease of exposition, step 1), where a warm-start is calculated, is described afterward.Step 2) is subdivided into two steps.Step 2.1) focuses on the first node and the first-time step/nodes of the second stage.This corresponds to problem CEP-A.The algorithm uses the current approximation of V 1 (x, x * , ω) to solve problem CEP-A (which approximates CEP) and obtain a candidate expansion plan x i , x * i .The algorithm proceeds with step 2.2), which aims at refining the approximation of V 1 (x, x * , ω) around x i , x * i .To do so, the algorithm performs a forward pass (step 2.2.1) and a backward pass (step 2.2.2).The forward pass proceeds forward in the number of second-stage nodes, solving the first-node subproblem and proceeding until arriving at the last node.The algorithm continues with the backward pass.This step computes supporting hyperplanes around the storages found during the forward pass and around the trial expansion plan x i , x * i .Starting from the last second-stage node, the subproblem is solved, using the dual multipliers for estimating a supporting hyperplane for the subproblem of the preceding node.The process is repeated until reaching the node associated with the first node of the second stage.The algorithm performs this procedure for all uncertainty realizations.Having described step 2), we can now describe the warm-start of step 1).Given an initial candidate expansion plan x 0 , x * 0 the objective is to provide an approximation of V 1 (x, x * , ω) around the given initial point.To do so, the warm-start performs step 2.2) throughout several iterations.These several iterations are performed in order to ensure that the storages found during the forward pass provide a reasonable approximation.
Remark 4.1: Note that both steps 2.1) and 2.2) can be solved efficiently.The former is a relatively small problem, which can be solved either by Bender's decomposition or subgradient schemes.In this work, we use Bender's decomposition.The latter step involves solving several small subproblems and thus does not pose a computational burden.As a consequence, this Algorithm 1 Second-stage relaxation algorithm.
(2.1) Solve CEP-A using the current approximation of V 1 .Thus, obtaining trial action x i , x * i .(2.2) for ω ∈ Ω: (2.2.1) Forward Pass.for k = 1, . . ., K: Solve the approximated problem 2) Using the dual multipliers compute a supporting hyperplane around x i , x * i , bv algorithmic approach is able to perform far more iterations than the subgradient algorithm.
Upper and lower bounds can be obtained as follows.An upper bound is computed by using the optimal values found during step 2.1) and during step 2.2.1).A lower bound can be obtained by the solution of step 2.1).By comparing upper and lower bounds, one can measure the optimality gap which can be used as a stopping criterion.The convergence of the decomposition algorithm is guaranteed due to the following lemma.
Lemma 4.2: Algorithm 1 converges in a finite number of iterations to CEP.
Proof: We prove that if no new cuts are added to V 1 (x, x * , ω) for ω ∈ Ω, then we are at the optimal of CEP.Suppose that no new cuts are obtained after iteration i and consider the expansion plan at that iteration x i , x i * .Note that, if no new cuts are added, then ω) (otherwise during the i + 1 iteration we would have found a new cut).Now let us assume that we have not converged, therefore there exists x, x * for which is an under approximation of CEP.Thus, its maximum possible value is the left-hand side of the previous inequality.As a consequence, This leads to a contradiction, therefore we conclude that x i , x * i is optimal.
Finally, let us show that this can be achieved in finite iterations.Following the same idea as in [48], where it is shown that SDDP-type [36] algorithms converge, we start from the subproblem of the last block.Note that the set of dual multipliers of the subproblem associated to corresponds to the vertices of the feasibility set of the dual problem, and it does not depend on x, x * , bv Consequently, there are finitely many supporting hyperplanes for the last subproblem, and thus after finitely many iterations we arrive at the supporting hyperplane description of V K (x, x * , bv t K−1 ,ω , v t K−1 ,ω , ω).We can now use such a description for the subproblem at block K − 1, and apply the same argument.Inductively, we can proceed backward in the number of blocks.

V. PARALLEL SCHEME
The present section aims at describing a parallel scheme for the algorithms presented in section IV.In this work, we have considered synchronous parallel implementations.
r Parallel subgradient algorithm: The parallelization strat- egy followed for the subgradient algorithm is presented graphically in Fig. 7.Each iteration proceeds as follows.
The master CPU provides an initial candidate expansion plan.The second-stage subproblems are distributed among the available CPUs.Each CPU solves its associated subproblem, the dual multipliers are collected and sent to the master CPU.At this point, the CPUs synchronize, i.e, a CPU stays idle until all other CPUs have finished their job.The master CPU uses the dual information in order to apply a projected subgradient and update the trial expansion plan.
r Parallel L-shaped algorithm: The parallelization follows a similar strategy as in the previous scheme.The difference is that during each iteration the master CPU calculates problem M.
r Parallel second-stage relaxation algorithm: The parallel scheme for our subgradient-relaxation algorithm follows a similar procedure.Step 2.1) is parallelized using the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
strategy described in Fig. 7. Step 2.2) is parallelized as follows.The uncertainty realizations are distributed among the available CPUs.Each CPU performs steps 2.2.1 and 2.2.2 of Algorithm 1 (see Fig. 6).At this point, the CPUs synchronize, so that the fastest CPU waits for the slowest CPU.

VI. CASE STUDY
The EVA problem aims at determining the expansion/retirement plans that will occur, and our study considers the 2025 target year.We have considered a total of 12 blocks per day.
Our algorithms are implemented in Julia v1.5 and JuMP v22.The chosen linear programming solver is Gurobi 9.The computational work is performed on the Lemaitre3 cluster of UCLouvain, which is hosted at the Consortium des Equipements de Calcul Intensif (CECI).The cluster consists of 80 compute nodes with two 12-core Intel SkyLake 5118 processors at 2.3 GHz and 95 GB of RAM (3970 MB/core), interconnected with an OmniPath network (OPA-56 Gbps).

A. Value of the Stochastic Solution
Due to the scale of the model, ENTSO-E has considered an approximate solution in ERAA 2021.This approximation proceeds by solving the so-called wait-and-see solution [18] of problem CEP.This leads to a candidate expansion plan x ω , x * ω , for each ω ∈ Ω.The average expansion plan x, x * is used as an approximation of the stochastic expansion plan.A natural question is whether this is a good approximation.To measure this, we can use well-established bounds.The following relationship holds [18]: Note that the right-hand side is the objective of CEP when using the sub-optimal average expansion plan x, x * .We refer to this as the average-W.S. solution.Note that these bounds do not require the calculation of the stochastic solution, and thus provide a reasonable way to measure if a stochastic solution is of interest for the problem.The wait-and-see solution has a cost of 5.0220e10 €, while the average wait-and-see solution has a cost of 5.2298e10 €.As one can observe, the relative difference between both quantities is approximately 4.13%, and thus the stochastic solution stands to improve the deterministic approximation by at most this quantity.This difference is of interest: adequacy studies such as ERAA aim at capturing adequacy metrics that involve curtailments.These curtailments, calculated as in (31), represent less than 2% of the total costs, and thus a solution whose total costs differ in 4% is clearly of interest.

B. Stochastic Solution
The previous subsection presented the relevance of a stochastic solution.The present subsection discusses how to obtain such a solution using the previously discussed algorithms.The algorithms have been run using 35 CPUs.In the case of the Fig. 8. Optimality gap evolution of the L-shaped method and our subgradientrelaxation algorithm.The x-axis is the elapsed time, while the y-axis is the relative difference between the upper and lower bounds.subgradient relaxation, we have decomposed the second stage into 92 blocks.
We begin by examining the convergence evolution of the L-shaped method.To do so, we examine the optimality gap evolution, which is presented in Fig. 8.As one can observe, the L-shaped method struggles to close the optimality gap.In fact, after more than a day of computation, the obtained gap is not of practical use.On the other hand, our subgradient-relaxation scheme is able to provide an optimality gap near 1% after about 4 hours of computation.
The proposed subgradient algorithm does not provide a lower bound estimate, therefore the calculation of an optimality gap is not presented.Fig. 9 presents the upper bound evolution of the subgradient algorithm and the subgradient-relaxation scheme.Recall that the subgradient-relaxation scheme uses an approximation of the economic dispatch, thus it does not correspond to the true value of using the expansion plan.Consequently, in order to have comparable upper bound values, after running the subgradient-relaxation algorithm we have evaluated the true cost of using the obtained expansion plan, this corresponds to the horizontal non-dashed green line.This quantity and the upper bound of the subgradient algorithm are comparable.The left panel of Fig. 9 presents the results when using 12 blocks per day, while the right panel presents the results when using 24 blocks per day.We note that both algorithmic schemes are able to converge.However, it can be observed that the subgradient-relaxation scheme converges considerably faster.In fact, for the 12-block, system after 30 hours of computation, the subgradient scheme has not attained the bound that the subgradient-relaxation scheme is able to find in just 4 hours.For the 24-block system, after almost 2 days of computation, the subgradient scheme fails to attain Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the bound that the subgradient-relaxation scheme finds in 15 hours.This indicates that the subgradient-relaxation scheme is better suited as the size of the problem increases.We note that the obtained run times can be used in practice.ENTSO-E's experience with a stochastic model, which consists of climate years and is solved as a large LP for a variety of target years, is that it requires longer run times.
Alternative approaches are considered for solving the problem, as a means of further comparison.This includes the progressive hedging algorithm [24] and Benders Decomposition With Multiple Master Problems (BDMM) [25].Unfortunately, these schemes are not suitable for our setting.On the one hand, each subproblem of progressive hedging is harder than Benders as (i) it includes a quadratic term in the objective, and (ii) includes the optimization of the first and second stage.This renders each iteration prohibitably expensive.On the other hand, BDMM increases the amount of second-stage Benders subproblems.In the case of the problem that we are interested in this is not tractable, as each Benders subproblem is large.
We highlight that the use of parallel computing has been crucial in order to obtain the solutions that are indicated here.In fact, as each CPU is handling one of the 35 climatic years, we expect a serial run to be approximately 35 times more time-consuming.

C. Solution Analysis
The present subsection aims at studying the differences between the stochastic solution obtained in subsection B and the deterministic solution obtained in subsection A. In terms of total costs, as presented in Fig. 10, we observe a difference of nearly 2.7% between the stochastic solution and the wait-and-see solution, a difference that supports the value of having computed a stochastic solution.Furthermore, we observe that the stochastic solution provides an improvement of nearly 1.5% with respect to the average W.S. solution, implemented by ENTSO-E in ERAA 2021.
This difference in total costs is examined in Fig. 11 by decomposing the total costs into operational costs, retired capacity savings, expansion plan costs, and curtailment costs.The upper panel presents the costs for both approaches, the stochastic solution, and the average W.S. solution, while the lower panel presents the relative difference between these solutions (a positive number indicates that the average W.S. solution has a higher value).Note that the stochastic solution tends to result in significantly fewer curtailments, which can be explained by the increased investments and fewer retirements, that is to say, we arrive at a more conservative solution.On the other hand, the lack of information regarding the future possible outcomes in the average expansion plan approach makes the solution prone to over-retiring and underestimating the required investments, which in turn implies that the power system is unable to satisfy the demand more effectively.

D. Adequacy Metrics
One objective of the ERAA study is to measure the ability of the system to maintain secure levels of supply.Two adequacy metrics are of particular interest to ENTSO-E.The first one is the loss of load expectation (LOLE), defined as LOLE = E ω [LOL ω ].Here, LOL ω is the number of hours during which demand is not served for the climatic year ω ∈ Ω.The second metric is the expected value of energy not served (EENS), which is defined as EENS = E ω [ENS ω ].Here, ENS ω is the number of curtailments obtained in climatic year ω ∈ Ω.We consider the 35 climatic years that are used for formulating the stochastic model, thereby aiming to compare its performance against the approximate solution that is obtained by averaging the expansion plans (which is the approach used by ENTSO-E in ERAA 2021).Therefore, there is no out-of-sample testing.We calculate the metrics of interest and present the results in Fig. 12.We consider both the LOLE and the EENS for the 4 main European regions: Central and Eastern Europe, Western Europe, Southern Europe, and Northern Europe.Fig. 12 effectively demonstrates how the use of a stochastic solution is able to provide consistent and significantly more accurate metrics, which can lead to a betterinformed adequacy assessment in Europe.This attribute is of particular interest to studies such as ERAA.

VII. CONCLUSION
In this work, we present a high-performance computing approach for obtaining a stochastic solution for the EVA of ERAA 2021.We specifically propose a subgradient algorithm implemented in parallel computing infrastructure, as well as a subgradient-relaxation approach that emerges as being practically attractive due to its capability to tackle large-scale problems efficiently.We further observe that the commonly used L-Shaped method is unable to provide a solution of practical use in a reasonable amount of time.Finally, we compare the stochastic solution against a deterministic approach implemented by ENTSO-E for the ERAA 2021 edition.A noticeable difference between both solutions in terms of commonly used adequacy metrics emerges, which highlights the practical value of the stochastic solution that we are able to compute.
Future ERAA studies aim at improving the EVA on two fronts.(i) On the one hand, a multi-year stochastic model is considered, that decides investments and retirements before the realization of each year.(ii) On the other hand, a more robust model is considered, which includes further sources of uncertainty, such as random outage patterns.These increased uncertainty motivates the study of scenario reduction techniques.
Furthermore, fuel price assumptions can have a significant impact on the results, and are also quite volatile, as the recent energy crisis suggests.The proposed approach can handle objective function uncertainty, thus allowing us to include these assumptions directly into the model.On the other hand, ENTSO-E's interaction with industry stakeholders points that this uncertainty can be considered as a sensitivity run of the model with different assumptions.This can provide insights on the impact without introducing further computational complexity to the model.

Fig. 1 .
Fig.1.ERAA methodology consists of 2 steps.The first step, the so-called economic viability assessment, is calculated using a set of uncertain climatic conditions.The second step computes adequacy indicators, which use the previously computed expansion plan over an uncertainty set that consists of climatic conditions as well as random outages.

Fig. 3 .
Fig. 3. Overview of climatic years of the four main European regions.The upper and lower bounds represent the 0.75 and 0.25 quantiles.The upper left panel presents wind production, the upper right is PV production, while the lower figure corresponds to hydrological inflows.

Fig. 4 .
Fig. 4. Decomposition of the two-stage stochastic problem.Each node represents an uncertainty realization and has associated with it an optimization problem that aims at minimizing the costs of the stage.

Fig. 5 .
Fig. 5. Chronological decomposition.The second stage is partitioned into several consecutive chronological blocks.

Fig. 6 .
Fig. 6.Step 2 of the second-stage relaxation algorithm.Step 2.1 solves the current approximation of problem CEP.Step 2.2 refines the approximation of problem CEP around the trial expansion plan found during step 2.1.

Fig. 7 .
Fig. 7. Parallel subgradient algorithm.The master CPU provides a trial expansion plan, which is sent to the subproblems.The subproblems are distributed among the CPUs and solved, the dual multipliers are sent to the master CPU.

Fig. 9 .
Fig. 9. Convergence evolution of the subgradient algorithm and the subgradient-relaxation scheme.The x-axis is the elapsed time while the y-axis is the cost.The yellow horizontal lines correspond to the wait-and-see solution and ENTSO-E's solution.The left panel corresponds to a setting with 12 blocks per day, while the right panel uses 24 blocks per day.

Fig. 11 .
Fig. 11.Cost breakdown.The upper panel presents the total costs, while the lower panel presents the relative difference between the quantities shown in the upper panel.