A two-stage stochastic MILP model for generation and transmission expansion planning with high shares of renewables

This paper is concerned with the generation and transmission expansion planning of large-scale energy systems with high penetration of renewable energy sources. Since expansion plans are usually provided for a long-term planning horizon, the system conditions are generally uncertain at the time the expansion plans are decided. In this work, we focus on the uncertainty of thermal power plants production costs, because of the important role they play in the generation and transmission expansion planning by affecting the merit order of thermal plants and the economic viability of renewable generation. To deal with this long-term uncertainty, we consider different scenarios and we define capacity expansion decisions using a two-stage stochastic programming model that aims at minimizing the sum of investment, decommissioning and fixed costs and the expected value of operational costs. To be computationally tractable most of the existing expansion planning models employ a low level of temporal and technical detail. However, this approach is no more an appropriate approximation for power systems analysis, since it does not allow to accurately study all the challenges related to integrating high shares of intermittent energy sources, underestimating the need for flexible resources and the expected costs. To provide more accurate expansion plans for power systems with large penetration of renewables, in our analysis, we consider a high level of temporal detail and we include unit commitment constraints on a plant-by-plant level into the expansion planning framework. To maintain the problem computationally tractable, we use representative days and we implement a multi-cut Benders decomposition algorithm, decomposing the original problem both by year and by scenario. Results obtained with our methodology in the Italian energy system under a 21-year planning horizon show how the proposed model can offer professional guidance and support in strategic decisions to the different actors involved in electricity transmission and generation.


Introduction
Generation and transmission expansion planning models aim at determining the least-cost investment schedule for constructing new power generation capacity, building new electrical interconnections and decommissioning thermal power plants. The definition of joint expansion plans is one of the most relevant problems in the field of power systems. Indeed, this kind of analysis provides a lot of useful information, allowing for instance to study the impact of some policy decisions and the possibility to achieve targets such as decarbonisation, integration of large shares of renewables or reduction of CO 2 emissions. Many different agents are interested in the joint generation and transmission expansion planning. Indeed, in all countries where the unbundling of the energy sector is still ongoing, this problem is addressed by vertically integrated monopoly utilities to establish a strategic master plan, to secure long-term supply, promoting affordable and reliable electricity and reducing power outages. Instead, in an unbundled energy environment, expansion planning analysis is performed by transmission system operators (e.g. TSO) for what is called anticipative planning, in which those planners determine the expansion plan that is optimal for the energy system as a whole to identify the best network reinforcement and to set incentives that could induce generation companies to invest in a socially efficient manner [1][2][3].
Since generation and transmission expansion plans are usually provided for a long-term planning horizon, the system conditions are generally uncertain at the time the expansion plans are decided. Different sources of uncertainty may affect planning decisions and must be considered in the decision-making process. These uncertainties can be classified into short-term and long-term uncertainties. Specifically, short-term uncertainties include the stochastic production from intermittent renewable energy sources and the demand variability throughout the hours of the day and the days of the week. These uncertainties are also known as random uncertainties, since they have an underlying probability distribution, which can be approximated from historical data. Instead, long-term uncertainty refers to long-term dynamics, including future values of investment costs, fossil fuel prices and policy constraints such as carbon prices. These uncertainties are also known as non-random uncertainties, since they are not usually probability-based: their probabilities are generally obtained based on expert judgement.
To accurately represent the uncertainty framework in the planning decisions, several approaches have been developed in the literature, including a fuzzy decision approach, chance-constrained models, Monte Carlo simulation, robust optimization, adaptation programming, and stochastic programming. In particular, given the relevant penetration of intermittent renewable power sources, many recent studies focus on the influence of these uncertainties on generation and transmission planning. For instance, in [4] an algorithm is developed for multi-objective optimization transmission expansion planning considering wind farm generation and combining Monte Carlo simulation and Point Estimation Method to investigate the effects of network uncertainties. Reference [5] proposes a chance-constrained formulation to tackle the uncertainties of load and wind turbine generators in transmission network expansion planning. Load and wind power generation are considered the main sources of uncertainty also in [6], where authors propose an efficient approach for probabilistic transmission expansion planning, dealt with by a Benders decomposition algorithm combined with a Monte Carlo method.
Other studies deal with long-term uncertainties. For instance, authors in [7] address the problem of transmission expansion planning under uncertainty in an electric energy system, considering future demand growth and the availability of generation facilities as main uncertainty sources. A robust optimization model is used to derive the investment decisions that minimize the system's total costs by anticipating the worst-case realization of the uncertain parameters within an uncertainty set. Instead, reference [8] proposes a robust generation and transmission expansion planning model, including flexible AC transmission systems (FACTS) devices and considering the uncertainty related to the annual net load duration curve. A robust optimization approach is also adopted in [9,10] to address generation expansion and transmission expansion, respectively. Specifically, reference [9] introduces a multiyear robust methodology for expansion planning, modelling the uncertainties associated with forecasted electricity load demand, as well as estimated investment and operation costs, through distribution-free bounded intervals producing polyhedral uncertainty sets. In [10] two optimization criteria for the transmission expansion planning problem under the robust optimization paradigm are studied, where maximum cost and maximum regret of the expansion plan overall uncertainties are minimized, respectively.
Adaptation programming represents another method to cope with uncertainty. As described in [11], adaptation programming designs a flexible system by minimizing the sum of investment and operational cost and system future adaptation cost to the conditions of other identified scenarios. Reference [12] further explores the adaptation programming method by applying this kind of model to a small system over a 40-years planning horizon, considering wind and solar investment cost, carbon taxes, demand and peak demand growth, fuel prices and transmission costs as uncertainty sources.
Among all the techniques developed to include uncertainty in the expansion planning framework, the most widely used is stochastic programming [13], a methodology introduced in the 1950s that uses a set of scenarios to model the future realization of the uncertain parameters in the considered planning horizon. In recent years, several studies in the field of stochastic programming have been carried out, leading to the development of two classes of methods: two-stage and multi-stage models. In a typical two-stage stochastic model, the investment decisions represent first-stage decisions, which are made before any uncertainty is revealed. Operational decisions are instead second-stage decisions, made after the realization of parameter values. For instance, in [14] a two-stage stochastic programming model for joint generation and transmission expansion is presented, considering as random events the demand, the equivalent availability of the generating plants and the transmission capacity factor of the transmission lines. Reference [15] presents a stochastic two-stage optimization model to evaluate interregional grid reinforcements in Great Britain. The same approach is also used in [16] to determine the type and quantity of power plants to be constructed in each year of an extended planning horizon, considering the uncertainty regarding future demand and fuel prices. Authors in [17] propose a two-stage stochastic generation expansion model, where the long-term wind power uncertainty is represented by a set of scenarios.
The two-stage approach can be extended to a multi-stage method, constructing models that are both more flexible and complex. As explained in [18], in multi-stage approaches expansion decisions are made at different stages, i.e., different points in time of the planning horizon. The expansion decisions at each stage depend on the scenario realization of the previous periods, but they do not depend on the future scenario realizations. Examples of multi-stage models are represented by [19][20][21]. For instance, in [19] a multi-stage multi-scale linear stochastic model is presented to optimize electricity generation, storage and transmission investments over a long planning horizon. Both long-term uncertainties, such as investment and fuel-cost changes and long-run demand-growth rates, and short-term uncertainties, such as hour-to-hour demand and renewable-availability uncertainty, are considered in this analysis and the progressive hedging algorithm is applied to decompose the original model by scenario. Reference [20] deals with wind power investments considering three major issues: the production variability and uncertainty of wind power facilities, the eventual future decline in wind power investment costs and the significant financial risk involved in such investment decisions. Recognizing the previous important issues, this paper proposes a risk-constrained multi-stage stochastic programming model to make optimal investment decisions on wind power facilities along a multi-stage horizon. Finally, in [21] authors study how uncertain future renewable penetration levels impact the electricity system and try to quantify effects for the Central European power market, by applying a multi-stage stochastic investment and dispatch model to analyse the effects on investment choices, electricity generation, and system costs. Although the multi-stage approach better represents long-term dynamics than the two-stage method, the complexity of the problem in multi-stage models is further increased. Finding the right balance between modelling accuracy and computation tractability remains an open research topic. Specifically, to limit computational costs most of the existing planning models employ a low level of temporal detail. However, as shown in [22], the low level of temporal detail can significantly affect the results, overestimating the renewable capacity and, thus, resulting not suited to accurately study all the challenges related to integrating high shares of intermittent energy sources. Moreover, most of the existing models for power generation and transmission expansion planning consider also a low level of technical detail, generally modelling thermal capacity expansion through continuous variables and, thus, without considering unit commitment constraints on a plant-by-plant level. However, as discussed in [23], this approach, although very common in the literature, is no more an appropriate approximation for expansion planning models. Indeed, due to the increasing penetration of intermittent renewable energy sources in power systems, the need is growing for flexible resources that could respond to the variability and uncertainty of stochastic generation. Ignoring unit commitment constraints leads to the impossibility of properly evaluating such flexibility and, consequently, to underestimate the required new generation capacity as well as the system costs. Table 1 reports a list of important features for expansion planning models, such as the modeling approach, the inclusion of investment decisions in new generation, storage and transmission facilities, the temporal detail in power system's operation evaluation, the inclusion of unit commitment constraints, the longterm policies considered in the expansion planning framework and the number of periods in which investment decisions can be made along the planning horizon. According to these characteristics, some relevant works in the literature and the model proposed in this paper are compared.
Specifically, to provide reliable expansion plans for large-scale energy systems with high shares of renewables, in this paper we propose a two-stage stochastic model for joint generation and transmission expansion planning, being decommissioning of thermal power plants and investments in new generation, transmission and storage facilities first-stage variables and operating decisions second-stage variables. The distinct feature of the proposed model is a detailed representation of the short-term operation. The choice of a two-stage model rather than a multistage approach is due to computational restrictions. Indeed, as can be noticed in Table 1, existing multi-stage models represent power systems short-term operation with a lower level of accuracy to the proposed model. For instance, no longterm storage facilities or unit commitment constraints are included in references [19,21], while the proposed model provides an accurate modeling of short-term dynamics by considering the unit commitment constraints on a plant-by-plant level to properly evaluate all the challenges related to integrating high shares of renewables. Moreover, to reduce the computational burden, existing multi-stage models consider only a limited number of periods in which investment decisions can be made along the planning horizon. Our approach differs from multi-stage models by allowing investment decisions in each year of the planning horizon.
Contributions of this paper are as follows: • High level of temporal detail by applying a clustering analysis on different time series (i.e., load, wind capacity factors, and solar capacity factors) to identify a small number of representative days. • High level of technical detail, including unit commitment constraints on a plant-by-plant level and providing an accurate representation of both shortterm and long-term storage. • Modelling of long-term dynamics, determining the investment schedule along the planning horizon considering multiple investment periods and considering long-term policy goals and environmental targets. • Modelling of long-term uncertainty, considering different scenarios for production costs of thermal power plants. • Application of Benders Decomposition, decomposing the original model both by year and by scenario to maintain the stochastic problem computationally tractable. • Application to a real-case study related to the Italian energy system, to evaluate the achievement of policy goals set by the European Commission. To the best of our knowledge, there are no other works with the same characteristics in the literature. The remainder of the paper is organized as follows. In Sect. 2 we describe the clustering analysis performed to identify the representative days. The mathematical formulation of the expansion planning model as a two-stage stochastic mixed-integer linear programming problem is introduced in Sect. 3. Section 4 describes the decomposition algorithm. In Sect. 5 we present a real case study concerning the expansion of the Italian energy system. Finally, conclusions are drawn in Sect. 6.

Clustering analysis
To provide a better representation of the short-term operation while maintaining the problem computationally tractable, some energy planning models use a small number of representative days instead of considering every hour of the planning horizon. Different approaches have been proposed to identify the representative days. Some authors use simple heuristics, selecting the days that contain the minimum demand, the maximum demand or the largest daily demand spread [24]. Other works combine heuristic approaches with the random selection of some additional days [25,26]. More advanced methods apply clustering algorithms to group days with similar load, wind production or solar production into clusters [27][28][29]. For each group either the cluster centroid or a specific historical day is then taken as the representative day. Finally, some works select representative days by minimizing the difference between the load duration curve and the one constructed by using the representative days [30,31]. In this paper, we identify representative days through a hybrid approach based on the iterative application of the k-medoids algorithm and on the selection of the most suited number of representative days to be used by considering load duration curves. We refer the reader to [32] for a full description of the proposed approach, hereafter summarized.
Specifically, this analysis aims to extract representative days from larger sets of days, consisting of multiple time series of electricity demand and renewables capacity factors. Therefore, a very important feature to be considered in the selection process is the correlation between different time series. In our approach, we apply the following procedure, designed to capture the spatial correlation between time series: Vectors creation We consider the first year of the planning horizon and for each day d we define a vector V d containing values of load, wind capacity factors and solar capacity factors for every zone and every hour of the day. Therefore, the dimension of each vector V d equals 72 times the number of zones, while the number of vectors is equal to 365.
Normalization Since wind and solar capacity factors are parameters in the range [0;1] , to properly evaluate the distance between vectors, hourly zonal load profiles are normalized, by subtracting the minimum load and dividing by the difference between maximum and minimum load. k-medoids algorithm application The objective of the k-medoids algorithm is to compute k clusters to minimize the deviation between observations V d and their representative V * c , which it is also called quantization error or intracluster distance: with c being the index for clusters, d the index for days and D c the set of all days d within cluster c . While for k-means parameters V * c are computed as the mean of all the points within the cluster c , for k-medoids the representative of each cluster is a vector belonging to that group. As discussed in [29], since using as representative days specific historic days rather than clusters' means usually provide better results, in our approach, we suggest applying k-medoids rather than k-means. Representative days weighting In our approach, representative days are weighted using cluster sizes, which are the number of observations grouped within a specific cluster. Selection of the number of representative days To define the best number of representative days to be used we perform the k-medoids algorithm for increasing values of the parameter k (starting with k = 2 ), comparing the obtained results in terms of load duration curves. Specifically, for each partition, we compute for each zone the mean absolute percentage error (MAPE) between the original zonal load duration curve and the one obtained by using the representative days and their associated weights. We then compute the system average MAPE (i.e., the average between the MAPEs of zonal load duration curves) and we select the minimum number of representative days that allows obtaining a system average MAPE lower than an input threshold. Application of annual load growth factors Once the representative days for the first year of the planning horizon are determined, the representative days of the following years are derived by applying annual growth factors to load profiles.
The identification of representative days allows reducing computational time while maintaining a high level of temporal detail, but it introduces a strong limitation in storage operation modelling. Indeed, as explained in [33], this approach allows providing a good representation of storage operation within a day, but since it does not preserve the chronology among representative days, any energy storage system with a cycle longer than 24 h cannot be modelled with the highest accuracy. To provide an accurate representation of long-term storage, in this paper we follow the same approach described in [29] applying constraints that guarantee some continuity between representative days. These constraints, along with the other equations that allow identifying the expansion plans, are described in the following section.
Investments in the new generation, transmission and storage facilities as well as decommissioning of thermal plants are defined using a mixed-integer linear programming model. This section presents the modelling assumptions as well as the notation and the mathematical formulation of the optimization model.

Modelling assumptions
The mathematical model that defines investment and decommissioning decisions is based on the system costs minimization formulation. Thus, adopting a centralized approach, the objective of the proposed model is to determine the expansion plan that is optimal for the whole energy system under study. Bidding strategies and competition between generators are beyond the scope of this work. In developing countries, where the unbundling of the energy sector is still ongoing, the centralized approach reflects the decision making framework, since expansion decisions are taken by vertically integrated monopoly utilities. Instead, in an unbundled energy environment, generation expansion pertains to profit-seeking agents, while transmission expansion pertains to the transmission system operator, which is a welfare-focused agent. In this case, the centralized approach reflects the anticipative planning analysis performed by transmission system operators to identify the best network reinforcement and to set incentives that could induce generation companies to invest in a socially efficient manner. As regards to investment decisions, in this paper we manage investments in new generating facilities differently. Specifically, while we use continuous variables to model wind and solar capacity expansion, investments in new thermal plants are managed through a set of candidate projects and binary variables describing whether the units are realized or not. Indeed, while it is possible to build wind and solar plants of any capacity, thermal units usually present a specified size, not allowing the modelling of thermal capacity expansion through continuous variables. Moreover, binary variables are needed to model some logical conditions between investment decisions, such as the existence of optional, mandatory, mutually exclusive, associate or upgrade projects. For the same reasons, binary variables model investments in new transmission lines. We introduce binary variables to represent also the building of new hydropower plants, usually characterized by predefined sizes, while investments in batteries capacity are modeled through continuous variables.
To represent the transmission network, we adopt a transportation model, imposing the fulfilment of nodal balance equations and transmission flow limits constraints, without incorporating voltage variables and the relationship of real power transfers with the bus angle difference and line impedance. Finally, the proposed optimization model is based on a perfect knowledge of investment costs, while demand is supposed to be inelastic.

Notation
To describe the two-stage mixed integer linear stochastic programming model for the power generation and transmission expansion problem, the following notation is introduced.

Y
Set of years Z Set of zones M Set of macro-areas K E Set of existing thermal power plants Set of existing thermal power plants to be mandatorily decommissioned Set of existing thermal power plants that may be optionally decommissioned K C1 ⊂ K C Set of candidate thermal power plants to be mandatorily constructed K C2 ⊂ K C Set of candidate thermal power plants that may be optionally constructed AK j ⊂ K C j-th group of associate candidate thermal power plants J AK Set of groups of associate candidate thermal power plants MEK j ⊂ K C j-th group of mutually exclusive candidate thermal power plants J MEK Set of groups of mutually exclusive candidate thermal power plants L E Set of existing transmission lines L C Set of candidate transmission lines L Set of transmission lines ( L = L E ∪ L C ) L C1 ⊂ L C Set of candidate transmission lines to be mandatorily constructed L C2 ⊂ L C Set of candidate transmission lines that may be optionally constructed AL j ⊂ L C j-th group of associate candidate transmission lines J AL Set of groups of associate candidate transmission lines  Year that contains representative day c cl(t) Representative day that contains hour t yr(t) Year that contains hour t Map ot,t Injective map of each hour ot to a representative hour t   Loss coefficient for energy stored by battery

Parameters
Upper bound on battery b charge : Power output of hydro reservoir h in hour t in scenario w E h,t,w : Energy level of hydro reservoir h in hour t in scenario w sl h,t,w : Spillage from hydro reservoir h in hour t in scenario w Ê h,ot,w : Energy level of hydro reservoir h in hour ot in scenario w E in b,t,w : Charge of battery b in hour t in scenario w E out b,t,w : Discharge of battery b in hour t in scenario w E b,t,w : Energy level of battery b in hour t in scenario w x l,t,w : Energy flow on transmission line l in hour t in scenario w ENP z,t,w : Energy not provided in zone z in hour t in scenario w OG z,t,w : Over-generation in zone z in hour t in scenario w RES z,t,w : Renewable generation in zone z in hour t in scenario w b. Binary variables k,t,w : 1 : thermal power plant k is ON in hour t in scenario w ; 0 : otherwise k,t,w : 1 : thermal power plant k is started up in hour t in scenario w ; 0 : otherwise k,t,w : 1 : thermal power plant k is shut down in hour t in scenario w ; 0 : otherwise

Mathematical model
The expansion planning model can be formulated as the following two-stage stochastic mixed-integer linear programming problem: The objective function (3.1) comprises the seven terms below: �� are the second-stage costs, i.e., the operating costs.
Specifically, item 7 above considers for each representative day the sum of production costs, start-up costs, hydro and batteries operational costs and penalties for energy not provided and over-generation. Production costs are supposed to be linear functions of the power output, being CM k,y,w the slopes of these linear relationships.
CM k,y,w = O&M k + HR k (Fuel.Price fuel(k),y,w + co 2fuel(k) CO2.Price y,w ) and optional projects, respectively. The associate project constraints (3.8) indicate that a group of projects is subject to a single investment decision, that is, either all projects in the group AK j are built, or none. Inequalities (3.9) are the mutually exclusive project's constraints and they ensure that only one unit (or none) of the projects in each group ( MEK j ) is built. Inequalities (3.10) and (3.11) impose lower and upper bounds on solar and wind installed capacity, respectively. Constraints Specifically, let k1 = UP(k) denote the new project that replaces the existing unit k when it starts operating: building project k1 within year y ( + k1,y = 1 ) implies the permanent offline status of unit k ( k,t,w = 0 ). Inequalities (3.26) state that the power output p k,t,w above the minimum power output P k is either bounded above by P k − P k , if unit k is online ( k,t,w = 1 ), or zero if unit k is offline ( k,t,w = 0 ). Constraints (3.27) enforce consistency between the binary variables k,t,w , k,t,w , k,t−1,w and k,t,w that represent start-up, shut down and status in adjacent hours, respectively, in all hours of the representative days, except the first hour of the day. Constraints (3.28) enforce consistency between k,t,w , k,t,w , k,c,w 0 and k,t,w , for the first hour of every representative day, where the parameter k,c,w 0 represents the status of unit k at the beginning of representative day c under scenario w . In our numerical experiments, the values assigned to k,c,w 0 are determined by applying a classification tree trained on historical data. The description of this procedure is outside the scope of this work and we refer the reader to [32] for a detailed description of how the values of parameters k,c,w 0 are computed. The minimum uptime constraints (3.29) impose that unit k can be started-up at most once in an interval of MUT k consecutive time periods. For each representative day, the minimum uptime constraints are enforced on the hours in the range from MUT k to the final (the 24th) hour of the day, to keep the representative days separate in the description of future scenarios. The minimum downtime constraints (3.30) work similarly to the shut-down of thermal power plants.
Constraints that all the energy charged and discharged since the previous check point plus the total energy at the previous check point are within bounds. This is possible because, as a result of the clustering procedure that determines the representative days, we know the representative hour associated with each hour of the planning horizon. This information is included in the model using the set Map ot,t . We refer the reader to [29] for a complete description of these long-term constraints. Equations (3.41) enforce the equality between energy levels of each reservoir h at the beginning and the end of the planning horizon.
Constraints (3.42) − (3.47) model the operation of batteries. Specifically, Inequalities (3.42) − (3.44) impose upper bounds to charge, discharge and stored energy and enforce consistency between the values of the first-stage and second-stage variables. Energy balances (3.45) apply to all hours but the first of the representative days and they ensure that the energy stored by battery b at the end of hour t equals the energy stored at the end of hour t − 1 (reduced by the loss coefficient b ≤ 1 ), plus the energy injected in b (reduced by the coefficient in b ≤ 1 ), minus the energy released from the battery (reduced by the coefficient out b ≥ 1 ). Equations (3.46) impose the energy balance for the first hour of each representative day, while constraints (3.47) state the equality for each battery b between energy levels at the beginning and the end of each representative day.
Inequalities (3.48) restrict the energy flows on the existing transmission lines. Constraints (3.49) impose lower and upper bounds to the power exchanges among zones and enforce consistency between the energy flows on candidate transmission lines and the binary variables related to investment decisions, not allowing energy flows on candidate lines which have not been built ( l,y = 0 ). The zonal balance Eqs. (3.50) impose equality between energy sources and use in every zone and every hour. Indeed, the left-hand side of these equations represents the hourly energy sources of zone z (given by thermal, solar and wind generation, incoming energy flows, hydro generation and energy released by batteries) and the right-hand sides describe the energy uses (represented by the load, outgoing energy flows, pumping power and energy absorbed by batteries). The variables ENP z,t,w and OG z,t,w allow detecting and evaluating problems in the simulated system that can cause a mismatch between supply and demand. Inequalities (3.51) ensure the fulfilment of zonal reserve requirements provided by thermal power plants (either existing plants not yet decommissioned or newly constructed plants).
Target constraints (3.52) − (3.55) model conditions required to promote sustainable development of energy systems. Specifically, inequalities (3.52) impose limits on thermal power generation employing fossil fuels whose availability could be limited in time. These constraints are imposed for each macro-area m , each year y , each fuel f , and each scenario w. In particular, they are computed by multiplying the daily consumption of fuel f under scenario w in all the zones belonging to macroarea m by the weight of each cluster, to take into account the different occurrences of representative days. Inequalities (3.53) impose limits for thermal energy production due to CO 2 emissions and they present a structure very similar to constraints (3.52), as also, in this case, the total daily emission in each representative day is multiplied by the cluster's weight. Equation (3.54) compute the zonal hourly renewable production, by considering solar, wind and hydroelectric generation. Constraints (3.55) control the renewables penetration, forcing the total renewable generation in macroarea m in year y to cover at least ratio m,y of the total yearly demand for electricity. Finally, constraints (3.56) − (3.68) define the optimization variables.

Decomposition algorithm
Given the long-term planning horizon and the high level of temporal and technical detail, the proposed two-stage stochastic programming model results computationally intractable even for a small number of scenarios. To obtain a solution, in this work, we apply a multi-cut Benders decomposition algorithm. Specifically, Benders decomposition is a method introduced in the 1960s [34] that allows solving a linear programming problem with complicating variables in a distributed manner at the cost of iterations [35]. In recent years, this algorithm has been widely applied to two-stage stochastic programming models [36][37][38]. Indeed, given their particular structure, two-stage stochastic programming models are suited for Benders decomposition application, being the first-stage variables the complicating variables: fixing the first-stage variables, the stochastic model decomposes into a set of independent and easy to solve subproblems. In the literature, there are several examples of power systems planning models solved through Benders decomposition. For instance, in [39] an enhanced Benders decomposition algorithm for two-stage stochastic linear problems is presented and applied to a large-scale dynamic generation and transmission expansion planning model for the European power system. In [23], authors implement a Benders decomposition algorithm to solve a two-stage stochastic generation expansion model, whose first stage determines the long-term expansion and short-term unit commitment decisions, while the second stage models the real-time operation. In [40] a Benders decomposition algorithm is used to solve a network-constrained AC unit commitment problem under uncertainty. Finally, references [41] and [42] apply Benders decomposition algorithm to solve large-scale transmission expansion planning problems.
In our model, the first-stage variables represent investment and decommissioning decisions and include − k,y , + k,y , wind z,y , sol z,y , h,y , cap b,y and l,y . If these variables are fixed, the original problem decomposes into a set of independent subproblems, one per year and scenario, each representing the operation in the second stage. Benders decomposition replaces the two-stage stochastic problem with an iterative collection of smaller problems. At each iteration, the so-called master problem is solved first to determine suitable values for the first-stage variables. Once the investment schedule is determined, the subproblems are solved. The number of these subproblems equals the number of years of the planning horizon times the number of scenarios. Finally, the dual information of the subproblems is sent to the master problem employing a cut to update the master problem solution. The next paragraphs provide the formulation of the master problem and the subproblems, as well as a more detailed description of the implemented algorithm.

Master problem
As previously mentioned, the master problem aims to provide values of the firststage variables by solving at each iteration j the following mixed-integer linear programming model (referred to as master problem (4) in the sequel) whose optimization variables are − k,y , − k,y , + k,y , + k,y , sol z,y ,wind z,y , h,y , h,y , cap b,y , l,y , l,y and w .
subject to The objective function (4.1) includes investment, decommissioning and fixed costs and the auxiliary variables w that lower-approximate the operation cost under scenario w . Therefore, the solution of this master problem represents a lower bound for the optimal solution of the original problem. At each iteration j , once the master problem is solved, the optimal values of the objective function and auxiliary variables w are stored in the vectors z l,y . These parameters will be used to build constraints (4.2) in subsequent iterations, along with the optimal value of subproblems objective function z (j) y,w (5.1) and dual variable vectors (j) y,w , obtained from fixing constraints (5.3)-(5.9) in the subproblems at iteration j , as described in the next paragraph. Specifically, the constraints (4.2), referred to as Benders optimality cuts, tighten the feasible region of the master problem over iterations. While in the original Benders decomposition algorithm, a single cut is generated at each iteration [34], in our approach we implement a multi-cut strategy, generating at each iteration one cut per scenario. As observed in [23,43] and [44], also in our application, the multi-cut Benders decomposition showed a faster convergence than the mono-cut algorithm. Moreover, it is worth mentioning that in our analysis, the master problem contains only optimality cuts, while Benders feasibility cuts are not included. Indeed, due to the second-stage variables ENP z,t,w and OG z,t,w that model energy not provided and overgeneration (i.e., energy in excess), respectively, the subproblems are always feasible. Lower bound constraints (4.3) on variables w avoid the master problem being unbounded in the first iteration. Constraints (4.4) control investment and decommissioning decisions as in the original problem, while constraints (4.5) define optimization variables.

Subproblems
At each iteration j , for given values of the first-stage variables l,y the subproblem associated with year y and scenario w (referred to as subproblem (5) in the sequel) is formulated as follows.

subject to
Objective function (5.1) minimizes the operating cost of the system in year y under scenario w . Constraints (5.2) include all operating constraints in the original problem and are imposed for each hour t belonging to the considered year y . Equations (5.3)-(5.9) fix the complicating variables to values determined by the master y,w , respectively. Both these parameters are needed to add Benders optimality cuts (4.2) to the master problem. Finally, constraints (5.10) and (5.11) define the optimization variables.
It is worth mentioning that constraints (5.10) replace the original definition (3.62) imposed on thermal commitment variables, i.e., in the subproblems the second-stage binary variables that describe thermal power plants activation patterns (i.e., k,t,w , k,t,w and k,t,w ) are relaxed to be continuous variables in the interval [0, 1] . Indeed, one requirement for Benders algorithm convergence is the convexity of subproblems. Thus, to obtain convex subproblems, binary unit commitment variables are relaxed to be continuous. In this way, the Benders algorithm will converge to a solution z * LP , which is optimal for the relaxed problem (with continuous and binary investment decisions and continuous operation decisions), but that not necessarily feasible for the original investment problem. For this reason, once the algorithm reaches convergence, the investment decisions are fixed and the subproblems are solved as MILP models (i.e., by considering binary unit commitment variables) so as to obtain a "quasi-optimal" solution for the original problem z MILP . The solution obtained with this procedure is therefore feasible, but not necessarily optimal (i.e., z * LP ≤ z * MILP ≤ z MILP ). However, in this application, as common in the literature when solving real-scale power systems, it is impossible to solve the expansion planning problem up to optimality. We consider a reasonable optimality gap tolerance of 0.1% (i.e., z MILP −z * LP z * LP ≤ 0.001 ). Empirical results show how the two solutions z * LP and z MILP are very close, being the relative distance in the current application 0.02% and, thus, much smaller than the reasonable optimality gap tolerance. Indeed, in our approach, we are modelling thermal unit commitment decisions by using the equations described in [45] and [46], which tighten the original problem's feasible region by reducing the distance between relaxed and integer solutions.
The solution of all subproblems allows computing the following upper bound to the optimal objective function value of the relaxed problem (with continuous and binary investment decisions and continuous operation decisions) at iteration j

Solution algorithm based on Benders decomposition
Given a small tolerance value ε to control convergence, the Benders decomposition works as follows: 1. Initialization Initialize the iteration counter, set j = 1 . Set z (j) y,w < , the optimal solution has been obtained, go to step 8. Otherwise, update the iteration counter, set j = j + 1 and go back to step 1.

Subproblems final integer solution
For each year y and each scenario w , solve subproblems (5) replacing constraints (5.10) with (3.62), i.e., considering mixedinteger linear problems. The solution obtained is now feasible for the original problem.

Numerical study
In this section, we describe the analysis conducted with the proposed model on the Italian energy system to assess the achievement at 2040 of policy goals set by the European Commission. A period of 21 years (planning horizon from 2020 to 2040) is simulated on the system divided into seven market zones. The simulation assumptions hereafter summarized and described are based on public information provided by ENTSO-E, ENTSOG, Terna, GME and EERA [47][48][49][50][51]. Input data and numerical results of the case study are described in the following subsections.

Data for Italian energy system
The market analysis considers the Italian electric power system divided into seven interconnected market zones: North (ITn), Central-North (ITcn), Central-South (ITcs), South (ITs), Calabria (ITcal), Sicily (ITsic) and Sardinia (ITsar). The neighbouring countries are modelled as four equivalent areas: Montenegro (ME), Greece (GR), Tunisia (TN) and one single zone named Europe (EU) that summarizes the power flow at the northern Italian border. Each of these equivalent areas is characterized by a set of equivalent power units whose bidding can model a dynamic Import/Export with the Italian power system. Figure 1 shows the existing interconnections as well as a set of candidate new interconnections among which the least cost generation and transmission expansion tool can choose. For modelling short-term operation, we have applied the procedure described in Sect. 2 by fixing at 5% the threshold for the system average MAPE in load duration curve approximation, obtaining five representative days for each year of the planning horizon. The value of 5% for the threshold has been chosen after several tests performed on the Italian scenario, which showed that this choice was a good balance between computational costs and approximation accuracy. By modifying the threshold, the number of representative days to be used changes as well, as shown in Table 2. Figure 2 illustrates the five representative days in year 2020 for the North zone, which has the highest electricity demand. As can be noted, each of the representative days is characterized by 24 values for load, solar capacity factors, and wind capacity factors and by a specific weight.
Representative days for the following years of the planning horizon have been obtained by applying an annual average demand growth of 1%, as shown in Fig. 3.
As far as renewable installed capacity growth is concerned, a lower bound of 50% of RES penetration (calculated as the ratio between renewable production and expected demand) has been imposed for year 2040. Based on this lower bound set by the Italian Energy Plan [48], the tool can decide how much to install and what is the best generation mix; however, it has to respect the minimum and maximum targets of photovoltaic and wind capacity set for every expansion plan year, according to the Italian Energy Plan.
Regarding the thermal fleet, based on data provided by the Italian TSO, the existing set of power plants has been implemented together with all related technical Fig. 1 Existing and candidate interconnections in the Italian power system information and decommissioning details. The achievement of the challenging RES penetration targets requires the availability of an appropriate reserve margin that can be provided by thermal generation, therefore a set of candidate Combined Cycle Power Plants and Open Cycle Power Plants (both fuelled by natural gas) has been considered to analyse how the system needs in terms of generation capacity and flexibility can be satisfied. Table 3 resumes the above-mentioned assumptions.
Regarding fuel consumption, we have considered a CO 2 emission cap of 70 Mtons for each year of the planning horizon, according to levels set by the European Commission, aimed at reducing the impact of the electricity sector on greenhouse gas emissions. In terms of the storage system, we have considered the possibility of adding new pumping units in the southern regions of Italy to the existing hydro pumping fleet and/ or the possibility of including in the system Lithium-ion batteries (cheaper but with less storage capacity) or Sodium-ion batteries (a little more expensive but with higher storage capacity). Table 4 resumes the technical data related to candidate storage projects.
The investment cost is assumed to decrease in the period 2020-2040, as stated in [50], with an exponentially decreasing trend as shown in Fig. 4.
For the economic factors, the values related to Investment Costs (IC), Fixed Costs and Decommissioning Costs are shown below, for both conventional thermal power plants and photovoltaic and wind power plants (see Table 5).
Both PV and wind investment costs are supposed to decrease along the planning horizon due to technology development. Specifically, while PV investment cost in the current analysis goes from 1000 €/kW in 2020 to 600 €/kW in 2030 to 450 €/kW in 2040, the wind investment cost is supposed to decrease from 1300 €/kW in 2020 to 900 €/kW in 2030 to 690 €/kW in 2040.
Fuel prices, together with CO 2 price, play an important role in the generation expansion plan optimization because they affect the merit order of thermal power plants and the economic viability of renewable generation. Moreover, fuel and CO 2 prices are quite volatile, depending on market fluctuations, the government's policy, and the political situation. For this reason, the stochastic analysis has been performed on prices, focusing on both the CO 2 and the fuel price. Indeed, we have used prices scenarios described in [47], which provides two scenarios for fuel prices (i.e., A and B scenario for fuel prices) and three different scenarios for emission costs (i.e., low, medium and high CO 2 prices). These scenarios are summarized in Tables 6 and 7. It is worth mentioning that codes Gas, Gasoil, and Coal in Table 7 represent the fuels used by Italian thermal plants, while code EUmix refers to the expected fuel cost in the foreign countries interconnected with the Italian market zones. Specifically, this cost reflects the expected variation of the generation mix in foreign countries, considering several factors such as increasing renewable penetrations and gas consumption, decreasing coal consumption and nuclear phase-out.

3
A two-stage stochastic MILP model for generation and…

Results and discussion
We solved the proposed model on a computer with two 2.  iteration we solved the master problem up to optimality. Figure 5 illustrates the evolution over iterations of the multi-cut Benders algorithm, which converges in 41 iterations. Indeed, at iteration 41 the relative distance between upper and lower bound equals 0.95 ⋅ 10 −4 , satisfying the predefined tolerance. The total time needed to solve the problem is 19,032 s, corresponding to 5 h 17 min and 12 s. Table 8 provides more information about computational times, specifying the size and the solution time for the master problem and the subproblems at the last iteration of Benders decomposition algorithm. As can be noticed, two different types of subproblems are considered, namely the base case and the updated subproblem. Indeed, language extension GUSS works as follows. First, the base case, i.e., the model instance related to the first subproblem, is considered. After solving the base case, the update data for each subproblem is applied to the model. Then, GUSS communicates the changes from the previous model instance to the solver. This procedure not only reduces the amount of data communicated to the solver, but also, in the case of an LP model, allows the solver to restart from an advanced basis and its factorization, dramatically reducing computational times. Indeed, as can be observed in Table 8, while the base case solution requires almost two minutes, each of the updated subproblems is   Moreover, it is worth mentioning that, although in this analysis we considered a small number of scenarios, thanks to Benders decomposition and language extension GUSS, the proposed model is scalable: a higher number of scenarios would not exponentially increase computational times. Specifically, Table 9 provides the solution time observed by considering the same tolerance for Benders convergence (i.e., = 10 −4 ) and a different number of scenarios, randomly built. As can be noticed, the stochastic model with 30 scenarios is solved in about 20 h, a result compatible with time requirements in investment studies.
As explained in the previous section, since in each iteration subproblems are defined as linear problems, second-stage decisions k,t,w , k,t,w and k,t,w , related to the commitment of thermal power plants, could be infeasible for the original problem. For this reason, once convergence is reached, subproblems are solved as mixed-integer linear problems, obtaining the final solution, characterized by total system costs higher than the convergence value of the Benders decomposition algorithm. However, thanks to the tight formulation of thermal unit commitment constraints, this difference is very small. Indeed, in this problem instance, the total system costs in the relaxed solution are 403,729 M€, while the integer solution is only 0.02% more expensive, with the objective function value being 403,810 M€.
The system expected costs for the whole expansion planning period are shown in Table 10. As can be observed, there is a remarkable difference between the first-stage and second-stage costs: while the sum of investment, decommissioning and fixed costs represents 15% of total costs, operating costs account for 85% of total costs. Specifically, the most relevant cost for the system is related to the production costs of thermoelectric power plants, which include O&M and fuel consumption, representing 99.6% of second-stage costs and 84.6% of total costs. On the contrary, start-up costs of thermoelectric power plants have a small impact, being 0.1% of total costs. While investment, decommissioning and fixed costs are independent of the scenario realization, operation costs, being second-stage costs, depend on the considered scenario for CO 2 and fuel prices. Specifically, Table 11 describes how system costs vary depending on stochastic prices.
As can be noticed, for each fuel price scenario A and B, the three CO 2 prices scenarios present similar values of start-up costs and operation costs of hydro plants, while batteries operational costs slightly differ between scenarios. As far as thermal production costs are concerned, they significantly differ in the six scenarios. As expected, since CO 2 and fuel prices affect the thermal plants marginal production costs CM k,y,w , the higher these parameters, the greater the production costs, which vary from 307,024 M€ in the fuel price scenario A with low CO 2 prices to 382,048 M€ in the fuel price scenario B with high CO 2 prices. Table 12 shows the additional capacity of wind and PV installed to reach the RES penetration target in 2040. The RES installed capacity consists of 59.2 GW of PV and 17.97 GW of wind power: this unbalance may be explained by the lower costs of the PV technology with respect to the wind technology. As far as interconnection projects are concerned, new national and international cross border lines must be implemented in years 2025, 2029 and 2040 to better exploit the stochastic renewable energy sources and compensate for the decommissioning of some Italian thermoelectric power plants. The selected interconnections are listed in Table 13.
Moreover, the tool couples the installed RES capacity with energy storage systems, installing throughout the planning period 5.9 GW of batteries and 3.98 GW of pumping units. A list that summarizes the installed capacity according to technology and zone is reported in Table 14. As can be noticed, as regards to batteries capacity, the model suggests installing both Lithium-Ion batteries and Sodium-Ion batteries in all market zones, diversifying capacity.
In the list of thermoelectric candidate projects, ten CCGT power plants have been selected as thermoelectric expansion capacity, starting operation in 2025. On the contrary, the decommissioning of one oil and three old CCGT power plants has been planned for 2027 in Sardinia, as a result of the forecasted increase of the natural gas price and of the high RES penetration in this zone. The new thermal power plants introduced in the system are located in the North, Central-South and South zones. These new thermal power plants ensure the availability of energy reserve margins. Figure 6 reports the expected energy generation in 2040 divided by energy source, for each Italian zone. As can be noticed, since the RES production (i.e., from wind, solar and hydro sources) is 218 TWh/year while the load is 404 TWh/ year, the target of reaching 50% of renewable penetration by 2040 has been fully achieved. Half of the total solar expansion capacity, as well as most of the thermoelectric expansion capacity, has been installed in the North, due to its high electricity demand, which is almost 60% of the Italian electricity demand. Even though generation exceeds demand in some market zones, in 2040 Italy will import 17 TWh from neighboring countries.  Figure 7 shows the installed capacity at the end of the planning period, grouped by source and zone. As can be noticed, the model suggests installing large shares of PV capacity in all Italian market zones and especially in the North, while the  wind expansion is mainly located in southern regions, which are characterized by the highest wind capacity factors.

Conclusions
In this paper, we have introduced a two-stage stochastic mixed-integer linear programming model for power generation and transmission expansion planning. The main characteristics of the proposed model are: • Identification of the least-cost investment schedule for investment in generating facilities, transmission lines, and storage systems and for decommissioning of existing thermal power plants. • High level of temporal detail to accurately capture short-term volatility of intermittent energy sources. • High level of technical detail, integrating into the expansion model the unit commitment constraints on a plant-by-plant level and accounting for both short-term and long-term storage to provide an accurate representation of the system's operation. • Modelling of long-term uncertainty, considering three scenarios for CO 2 prices combined with two scenarios for fuel prices, therefore considering six scenarios for production costs of thermal units. • Consideration of policy goals such as CO 2 emission limits and renewable penetration targets in expansion planning.
To be computationally tractable, the optimization model has been integrated with a clustering analysis to select representative days considering the correlation between different time series (i.e., daily zonal electricity demands and renewable capacity factors). To further reduce computational costs, a multi-cut Benders decomposition algorithm has been implemented, decomposing the original twostage stochastic programming model by year and by scenario.
The proposed model presents many applications in power systems analysis, allowing in general to address what-if questions. For example, in this paper, we have applied the proposed model to evaluate the possibility for the Italian energy system to achieve by 2040 a 50% renewable penetration target and a sustainable reduction of CO 2 emissions, according to the goals set by the European Commission. Empirical results show how solar PV technology could play a key role in achieving these Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.