Aﬃne Decision Rule Approximation to Address Demand Response Uncertainty in Smart Grids’ Capacity Planning

,


Introduction
Generation expansion planning (GEP) is a classical energy problem that aims at determining the required power capacity to satisfy demand over a long-term horizon at minimum cost, while satisfying economic, environmental and technical constraints (see Koltsaklis and Dagoumas 2018, for a recent review). For this, one must consider existing and future electricity generation technologies to determine an optimal investment and retirement plan for the power sector. Computationally speaking, in a deterministic setting, the GEP problem is a large-scale optimization problem that typically accounts for more than hundreds of thousands of decision variables and constraints in order to identify realistic solutions. Motivated by advances in the structure of smart grids, recent versions of GEP problems have become even more complex as they attempt to model the notion of demand response (DR) (see e.g. Babonneau et al. 2017, Lohmann andRebennack 2017). DR can be defined as the shifting of demand by consumers, in response to supply-side incentives offered at times of high wholesale market prices or when system reliability is jeopardized (US Department of Energy 2006). Development of sophisticated data-gathering devices such as smart meters and voltage sensors, as well as increasing penetration of flexible loads such as electric vehicles (EV) and heating pumps, provide the technical possibility of implementing DR programs through utilities. In the wholesale capacity market, the participants bid on providing DR resources as an alternative for expensive generation and transmission expansion strategies. DR resources can drastically reduce marginal production costs, as they are substitutes for technologies that help meet peak loads. However, including DR in the GEP problem is challenging. Indeed, DR is a price-sensitive resource, whose availability is not fully known, especially in the capacity expansion planning phase. As the planning horizon in a GEP problem is typically several decades, any planning or prediction of the availability of DR resources will be affected by two sources of uncertainties: i) errors in the prediction of total demand; and ii) changes in the response behavior. Failure to properly consider DR uncertainty (DRU) is likely to lead to situations in which there are insufficient capacity resources to satisfy realized demands.
Robust Optimization (RO) and Stochastic Programming (SP) are the two main methodologies to address uncertainty. The underlying assumption in SP is that the probability distribution of an uncertain parameter is known (Birge and Louveaux 2011). Conversely, RO does not assume information about the distribution of an uncertain parameter, but instead assumes that uncertain parameters lie in a user-defined uncertainty set. RO aims at finding solutions that are immunized to all perturbations of the uncertain parameters in the uncertainty set (Bertsimas et al. 2011). In multi-stage versions of SP and RO, it is assumed that decisions are taken at different points of time. The simplest form, known respectively as two-stage SP or Adjustable RO (ARO), considers here-and-now decisions that need to be made before having any information about the uncertain parameters and wait-and-see decisions that can be adjusted according to the realized uncertain parameters at a later time. It is generally known that ARO (Ben-Tal et al. 2004) provides less conservative optimal solutions than RO, since it has more flexibility to adjust the decisions with respect to the uncertain parameters. However, the better performance comes at the price of problem intractability. To circumvent this issue, Ben-Tal et al. (2004) propose using affine decision rules for the delayed decisions, a technique also referred as the affinely adjustable robust counterpart (AARC). Following Ben-Tal et al. (2004), other types of decision rules have been proposed for ARO problems Georghiou 2018, Georghiou et al. 2020). For a review of the ARO approach, we refer the reader to Yanıkoglu et al. (2019).
A number of different multi-stage robust or stochastic approaches have been previously used in the literature to address uncertainties in GEP problems (Dehghan et al. 2014, Mejía-Giraldo and McCalley 2014, Domínguez et al. 2016, Amjady et al. 2018, Baringo and Baringo 2018, Han et al. 2018, Zou et al. 2018. Because of the inherent computational challenges that emerge with these models, there should be a compromise between the level of details and the size of the problem. For example, Han et al. (2018) propose a two-stage stochastic program that only considers a single-period planning horizon, while the multi-stage model proposed in Domínguez et al. (2016) is only implemented in a case study involving four periods. It therefore appears that tackling uncertainty in longterm industrial-size GEP problems is out of reach for currently available methods. Improving the computational efficiency of algorithms for solving such problems is also crucial in helping to understand better the interactions between demand and supply side, as such studies can require iterating through the resolution of a number of GEP problems to identify market equilibria (see Babonneau et al. 2016Babonneau et al. , 2020, for examples of this approach).
In this paper, we propose a numerical method for addressing DRU on a realistic large-scale GEP problem. To do so, we use the Energy-Technology-Environment Model (ETEM) proposed in Babonneau et al. (2017) to model the market mechanism that matches flexible load and supply, where we introduce for the first time DRU as an implementation error of the DR decision variable. We formulate the problem as a robust multi-period linear program. In a first stage, the policy-maker decides how to invest in the capacity of each generation technology and plans for a target DR for each future time period. Then, periodically, depending on the actual level of contribution in DR programs for this period, the planner decides on optimal energy procurement based on available resources. We evaluate the performance of the proposed robust ETEM model on a real-world case study based on the energy system of the "Arc Lémanique" region in Switzerland. To solve the problem, we derive a Robust Multi-Period Conservative Approximation (RMPCA) of the problem, and develop a Benders decomposition algorithm (inspired from Ardestani-Jaafari and Delage 2018) to solve it. In order to focus on the effects of DRU on the design of smart electricity networks, we have limited the model in this paper to only consider the uncertainty of DR. However, the proposed methodology could be extended to accommodate supply-side uncertainties, such as availability of intermittent resources, or investment costs. We however leave this as the subject of future work as discussed in Section 7.
Overall, the contributions of this paper can be summarized as follows: 1. We introduce DRU in a GEP problem, addressed using the ETEM energy model, and propose a robust multi-period linear program in which the procurement decisions are adjustable to the actual DR. The approach that we propose is flexible in the sense that, with minor modifications, it could accommodate other sources of uncertainties in GEP problems such as uncertainty of intermittent energy resources and investment costs.
2. We propose, for the first time, to seek Pareto Robustly Optimal (PRO) solutions (see Iancu and Trichakis 2014) of the master problem in a Benders decomposition algorithm in order to accelerate convergence. When combined with the valid inequalities proposed in Ardestani-Jaafari and Delage (2020), solution time is improved by 40% on average in our randomly generated instances. Note that this idea differs from the idea of Pareto optimal cut, introduced by Magnanti and Wong (1981), in the way that the latter seeks a Pareto optimal solution among multiple optimal solutions of the sub-problems. In contrast, we seek a PRO solution in the master problem, i.e., a solution that is guaranteed not to be dominated by any other optimal solutions in terms of its constraint margin profile.
3. Compared to previous work in the literature, we solve an adjustable robust multi-commodity GEP problem whose size, in deterministic form, is two orders of magnitudes larger than the largest instance addressed in previous case studies. Moreover, our case study confirms in a simulation that adjustable robust policies can significantly improve (i.e., around 33%) the expected total cost of the energy system compared to the solution of the deterministic version when accounting for electricity shortage penalties. Furthermore, it shows that the flexibility of energy procurement decisions can be responsible for a reduction of 9 billion Swiss francs (CHF) in expected total expansion planning costs.
The remainder of the paper is organized as follows. In the next section, we present an overview of different approaches to model DRU, and summarize studies that use a multi-stage approach to model uncertainties in GEP problems. In Section 3, we introduce briefly the ETEM model. In Section 4, we present how we have modeled DRU in ETEM and derive the multi-period robust optimization problem as well as its conservative approximation. In Section 5, the Benders decomposition is detailed. Section 6 presents the case study and provides numerical results. Finally, Section 7 provides concluding remarks.

Literature Review
In this section, we first review papers that use a multi-period approach to model a robust or stochastic GEP problem. We specifically summarize their solution method and the size of the instances they solve. Usually, the size of a GEP problem is dependent on the number of time periods, load duration curve (LDC) steps, technologies, and commodities. Hence, we will summarize the size of studied instances with the product of these quantities, referred as the Deterministic Size Indicator (DSI) index. In a second part, we review recently proposed strategies to model DRU in different energy-related problems. Finally, we review the literature on the ETEM model. Bloom (1983) and Bloom et al. (1984) are among the first papers that propose a two-stage stochastic programming formulation for GEP problems with uncertain demand and supply. They propose a general Benders decomposition, where a master problem accounts for investment decisions, and a set of subproblems represent the annual operation cost and reliability level of the installed capacities while enforcing probabilistic reliability constraints. Han et al. (2018) develop a two-stage stochastic program to model the uncertainty of load demand and wind output in a GEP problem. Using affine decision rules, they reformulate the problem as a deterministic second-order cone optimization problem. They solve an instance with 1 period, 5 LDC steps, and 32 technologies to test their algorithm (DSI=160). Mejía-Giraldo and McCalley (2014) propose an adjustable robust optimization framework to address the uncertainty of fuel price, demand, and transmission capacities in a GEP problem. In this setting, investment decisions are set as an affine function of fuel price, whereas voltage angle (decision) is parameterized as an affine function of demand and transmission capacities. They reformulate the problem as an LP and test it on a simplified version of the US power system with 20 periods, 3 LDC steps, and 13 technologies (DSI=780). Domínguez et al. (2016) use linear decision rules to reformulate a multi-stage GEP problem and cast it as a tractable LP. They solve a GEP problem with 4 decision periods, 18 technologies and 100 LDC steps (DSI=7200). Zou et al. (2018) propose a partially adaptive multi-stage stochastic GEP problem. In this setting, the capacity expansion plan is adaptive to the uncertain parameter up to a certain point. Fuel price and demand are two sources of uncertainties in this study. They test their algorithm on a case study with 10 periods, 3 LDC steps, and 6 technologies (DSI=180). Dehghan et al. (2014) develop a two-stage robust optimization problem in which the load demand and investment cost are uncertain. They solve their problem, using a cutting plane method. The case study is a 10-period planning of a network with 120 technologies and 3 LDC steps (DSI=3600). None of the above studies model neither DR nor DRU. In addition, all of them are a single-commodity GEP problem, meaning that they only model the electricity network and account for an aggregated electricity demand over the planning horizon. By contrast, in a multi-commodity energy model such as ETEM, one also models the interactions of the electricity sector with other sectors of the energy system (e.g., heat, transportation sectors), as well as the demands for final energy services (such as heat and lighting). This increases the size of the deterministic problem. In this paper, we solve a case study with 10 decision-making periods, 12 LDC steps, 142 technologies, and 57 energy commodities (DSI=971,280).
DR is also modeled in other power system planning problems such as unit commitment, transmission or distribution expansion planning, and wholesale electricity market problems. In general, there are two strategies one can follow to implement DR programs: dispatchable and non-dispatchable. In the former, the utility directly controls the timing of some loads of the customers who voluntarily participate in the program. More precisely, the utility cuts down the consumption during peak periods to ensure system reliability, and in return, remunerates participants with annual payments. In non-dispatchable strategies, the utility sends to customers a price signal with the purpose of flattening the demand curve. Automatic price-responsive devices adjust the timing of the consumption. While developing an optimal consumption schedule is the focus of some research , another stream studies the conditions under which the dynamic price signals improve social welfare (see, e.g., Adelman andUçkun 2019, Venizelou et al. 2018). Elaboration on the different strategies to implement DR and an overview on the success of these programs in the US electricity market is presented in Shariatzadeh et al. (2015).
The effect of DR on the wholesale and capacity market is another stream of research. Vatani et al. (2017) introduce DR to a capacity market problem in which the cost-minimizer capacity planner trades off between new generation and transmission capacity expansion and DR expenses. Arasteh et al. (2015) model a trade-off between DR expenses and distribution expansion costs. Zhang and Zhang (2019) model DR as an alternative procurement strategy for an electricity retailer. Namely, the electricity retailer must procure energy in a day-ahead wholesale market. Given that the prices in this market are highly uncertain, DR is introduced as a leverage to manage the retailer's risk. In this approach, the retailer solves its problem to determine a procurement strategy as well as DR incentive prices. Then, consumers respond to the signaled price in order to minimize their consumption cost. Babonneau et al. (2016) model DR in a capacity expansion planning problem considering distribution constraints. They develop a game framework between the retailer and many small consumers. The retailer solves a minimization model to determine his desired DR and capacity expansion plan. Then, using marginal cost of production, which corresponds to dual variables, he signals consumers to adjust their consumption. They argue that, if the price function, with which the retailer signals consumers, is convex with respect to both demand and required level of reserved capacity, their approach can be cast as a linear programming problem. Finally, Lohmann and Rebennack (2017) model DR in a long-term capacity expansion planning problem.
In all the above research, DR was mostly modeled in a deterministic way. In other words, when DR was modeled using a price elastic demand curve, no uncertainty affected this curve. Alternatively, when it was modeled as a decision variable, no implementation error was affecting optimal DR decisions. However, due to a variety of reasons, including demand prediction errors and changes in response behavior, in practice DR necessarily ends up different from what the supply side expected. According to Chatterjee et al. (2018), a lack of DR forecasting and estimation tools is one of the most important barriers to wider adoption of DR. On this regard, an active stream of literature has proposed a number of statistical learning methods for predicting demand response associated to storage-like flexible loads such as residential and industrial air conditioning (see, e.g., Dyson et al. 2014, Qi et al. 2020, and a stand-alone micro-grid with the presence of photovoltaic and wind generation units (see Amrollahi and Bathaee 2017). Another stream of research investigates the influence of DRU on its integration to the network. Li et al. (2015) model incentive-based DR as an alternative for transmission upgrade investments and consider consumer's bid for load reduction uncertainty. They develop a stochastic programming model and use Monte Carlo simulation and Benders decomposition for resolution. Ströhle and Flath (2016) model DRU in a local online market to match flexible load and uncertain electricity supply. Gärttner et al. (2018) formulate DRU in a demand aggregator's problem where the scheduling of demand and the optimal dispatch are optimized at the same time. He et al. (2019) develop a two-stage distributionally robust problem to model DRU in a distribution network expansion problem. Prices are modeled as a decision variable to optimally shape the desired demand response. However, the price elasticity is considered to be a random parameter. Moreover, distributed generation resources are also affected by uncertainty. They propose a column and constraint generation approach to solve the resulting problem. Zhao et al. (2013) model the uncertainty of DR and wind generation in a unit commitment problem. DR is modeled through a price-elastic demand curve with price uncertainty. They propose a three-stage robust optimization formulation. First stage decisions consider unit commitment planning, i.e., turning on and off the generators. Second stage decisions include the dispatch of electricity and are implemented after the actual wind power output is realized. Finally, in a third stage, the demand response is observed. Moreover, they propose a Benders decomposition to solve the three-stage problem. An overview of the different techniques to model ancillary services in unit commitment problems is presented in Knueven et al. (2020). Asensio et al. (2018) model the DRU through a price-sensitive demand curve in a distribution network expansion problem and use a scenario-tree approach to solve the resulting stochastic programming model. Roveto et al. (2020) develop a data-driven distributionally robust optimization model of the demand response auction, where the risk of the retailer is controlled using value-at-risk and conditional valueat-risk measures. They model DRU by considering DR bid uncertainty. In other words, the aggregator receives bids from consumers to reduce their consumption at a specific time. However, consumers might not actually reduce their consumption to this full amount. The aggregator aims to ensure a certain level of demand response procurement at minimum cost. Huang et al. (2020) develop a multi-stage stochastic programming model to optimally schedule power generation assets when the system operator has access to the ancillary service market. In this paper, we model DRU in a multi-commodity energy planning model, called ETEM (Babonneau et al. 2017), to address a GEP problem. In ETEM, DR is modeled through a decision variable which is optimized by the supply-side. It is assumed that the supply-side can shape consumer's response behavior through different long-term DR programs. For the first time, we model DRU as an implementation error of the DR decision in ETEM, and develop an adaptive robust multi-period conservative approximation formulation of the problem.
ETEM first was formulated by Babonneau et al. (2012). Although the developed formulation is deterministic, they also provide a stochastic programming model with different scenarios for the electricity import cost. As an extension to the primary ETEM model, Babonneau et al. (2016) propose a linear approximation of distribution system and power flow constraints in ETEM. The new formulation, called ETEM-SG, accounts for provision of reserved capacities provided by demand response in smart grids. Babonneau et al. (2017) and Babonneau and Haurie (2019) propose robust versions of ETEM to address the uncertainty of investment costs and of availability of technologies and transmission capacity of electricity transmission lines. In both papers, a static robust version of ETEM is developed and a relaxed version of the model is solved to reduce conservatism of the solution. Babonneau et al. (2020) is the first version of ETEM that employs robust optimization to handle uncertain demand response. Specifically, it proposes a robust ETEM model coupled with a mean-field game equilibrium model that represents the charging behavior of electric vehicle owners. The paper models the uncertainty of DR by introducing in ETEM two upper and lower bounding constraints on the sum of DR over all time periods. The bounding constraints are calibrated based on the confidence intervals of the random demand which is simulated using a mean-field game model. We note that Babonneau et al.'s way of accounting for DRU in ETEM is fundamentally different from ours. First, they again only involve static decisions and focus on a single-period investment problem. More importantly, a careful analysis of their approach reveals that it effectively offers no protection against the DR deviations in our setting. A more detailed demonstration of this weakness is presented in Appendix A. Our paper should therefore be seen as presenting a significant extension of ETEM that identifies adaptive long term strategies to immunizes for the first time the energy system against DR deviations.

The ETEM Energy Model
As proposed in Babonneau et al. (2012), ETEM is a long-term, bottom-up energy model cast as a linear programming problem. It is a member of the MARKAL-TIMES family of models, and represents the entire energy sector from primary resources to useful energy demands. ETEM is a demand-driven model in which the objective is to provide energy services at minimum cost. The planning horizon is typically more than 50 years and the model gives insights on both long-term strategic and short-term operational decisions. While the long-term decisions include capacity expansion and strategic targeting for a desired demand response, short-term decisions consist of energy procurement planning, which includes optimal generation, import, export and regional transmissions. ETEM can also account for technical, economic and market constraints (see Babonneau et al. 2017).
ETEM is based on the concept of a reference energy system (see Figure 1 for the depiction of the energy system of the Arc Lémanique region in Switzerland considered in our case study). The model considers different energy commodities (captured by vertical lines in the figure) namely "resources" (i.e., a mix of primary energy resources and some imported secondary energy forms), selected secondary forms (i.e., electricity and heat), and useful energy demands. More precisely, resources correspond to coal, oil products, gas, hydro, intermittent renewables (wind and solar), as well as other sources such as municipal solid waste and wood. Useful energy demands belong to four categories: industry, residential electricity, residential heat, and transportation. The boxes correspond to the different energy technologies, from generators to technologies providing energy services to end-users. In our implementation, we also add an energy resource (denoted "ExR") with a high cost that represents the cost of shortage in the system. A dummy technology (denoted "EXT") links this expensive resource to the useful demands (see Assumption 1 and following discussion for more details).
Figure 1: Arc Lémanique reference energy system, where Int-Res stands for intermittent resources (solar and wind), N-Gas for natural gas, Oil-fuel for processed oil products, Other for hydrogen and additional sources of energy (such as solid waste and wood and geothermal), ELC for electricity, RES for residential, PP for power plant, CHP for combined heat and power plant, IMP for electricity import, Ind-Machinery for industrial machinery, ExR for a dummy expensive resource added to avoid in-feasibility, ExT for a dummy expensive technology used to avoid infeasibility. Other-PP includes geothermal, fuel-cell and municipal waste power plants.
The planning horizon in ETEM spans over several decades to simulate the long investment cycles typical of the energy sector. We denote by t ∈ T := {1, . . . , |T|} the index representing the decision-making periods. In addition, to capture short-term operational patterns of the demand load, each period is divided into smaller time-slices s ∈ S, representing load duration curve steps. In ETEM, 12 time-slices are assumed (S = {1, · · · , 12}), which are partitioned among a set J of 3 typical seasons (winter, summer, and intermediate), and each day of these seasons is divided into 4 time slots (morning, midday, midnight, and night), as shown in Figure 2. Before describing ETEM's mathematical formulation, it is worth introducing our nomenclature. Let C be the set of energy commodities and let P be the set of energy technologies. The set of indices F identifies the different input-output energy flows associated to each technology. For instance, in Figure 1, Combined Heat and Power (CHP) uses natural gas, solid waste and wood as input to generate both electricity and heat. Finally, let L be a set of buses in different geographical zones. A full nomenclature is presented below.
Nomenclature t ∈ T Index for time period s ∈ S Index for time-slices p ∈ P Index for technologies c ∈ C Index for energy commodities c s ∈ CS Index for energy storage f ∈ F Index for energy flows l ∈ L Index for buses (geographical zones) j ∈ J Index for seasons i ∈ I Index for period-seasons (t, j) The objective function (1a) minimizes a discounted sum of all costs of the system over all regions (l ∈ L) and time periods (t ∈ T), where parameters α t,p , π t,p , ν t,p , λ t,s,c , λ t,s,c , and λ t,s,l,l ,c are unit costs of investment, fixed, variable, import, export, and regional transmission respectively. Constraint (1b) is a commodity balance constraint. It ensures that the regional procurement of each energy commodity c is greater, or equal than its total consumption at each period t and time-slice s. Specifically, the left-hand side of this constraint consists of i) total production of commodity c in region l by all technologies producing it (P P c ) ii) import of commodity c, and, iii) net transmission of commodity c into region l. It is worth mentioning that in ETEM, import (and similarly export) of energy refers to amount of input (output) energy to (from) the energy system from outside of it, but energy transmission refers to amount of energy which is produced in one region and transmitted to another region inside energy system. Parameter η c is the network efficiency with respect to commodity c, e.g., the efficiency of electricity transmission lines. On the right-hand side, the consumption is equal to the total consumption by technologies consuming the commodity, i.e., P C c , added to the amount of commodity that is exported. Constraint (1c) introduces a safety margin in procurement of commodity c ∈ C G , mostly electricity, during peak time-slices s ∈ S G , to protect against random events not explicitly represented in the model. Parameter t,s,c ∈ [0, 1] represents the fraction of reserved capacity needed to ensure covering the peak load. The left-hand side of this constraint models the maximum amount of commodity c ∈ C G that can be procured in period t and time-slice s. This amount is equal to the sum of i) the maximum production capacity of commodity c by technologies that produce it as their main output (P c ), ii) the production of commodity c by technologies that produce c as their by-product of their main activity, and iii) the import of commodity c. The left-hand side is the total consumption similar to constraint (1b). Parameter θ c p is the proportion of technology production that can be used during the peak period. Constraint (1d) models the balance of energy storage between consecutive time-slices. Amount of storage at time-slice c can be consumed at subsequent time-slice s . The notion of subsequency of time-slices is presented in Figure 2. Constraints (1e) -(1g) are the key constraints that model the use of demand response. Parameter Θ t,l,c is the total demand of service c ∈ C D in the period t and zone l. Variable V t,s,c is the demand response which optimally distributes the total demand of period t into all time-slices s inside period t. Constraint (1f) limits the demand response to vary between an interval around the nominal value, i.e., υ t,s,c . In addition, the sum of the demand response must be equal to the sum of the nominal values in each season, according to constraint (1g). Constraint (1h) limits the maximum production or consumption of each technology to the available capacity of that technology. The parameter β t,s,p is the capacity factor. The capacity factor of technology p is defined as the average power generated divided by the rated peak power. This simply represents the fraction of the total capacity which is available at each time-slice. In addition, since renewable generation is usually given priority in dispatch over conventional forms of generation (see for example Knueven et al. 2020, Lohmann andRebennack 2017), their production is imposed to take its maximum production capacity in constraint (1i). Constraint (1j) models the efficiency of technology p. In ETEM, the efficiency of a technology is modeled through the annual parameter η t f,f that links each output flow f ∈ FO p to an input flow f ∈ FI p . The constraint ensures that the sum of all output energy commodity c which is linked to output flow f , i.e., C f , is equal to a fraction of the sum of all input energy commodities c which is linked to input flow f , i.e., C f . We note that constraints (1k) and (1l) impose upper and lower bounds for seasonal production of technology p in region l, and seasonal import of energy commodity c in region l. Moreover, in constraint (1m) the seasonal procurement of commodity c, over all regions in L, is imposed an upper bound. This departs slightly from the original ETEM formulation, which only imposed annual bounds but is reasonable and enables the decomposition scheme presented in Section 5. We will discuss later that seasonally independent procurement decisions enable us to decompose the problem into several independent sub-problems where one decides on the optimal procurement given a capacity plan and demand response. This structure benefits the efficiency of the solution method. Finally, space X and Y represent other operational, technical, and economical constraints that define a desirable space for capacity and procurement decisions. Moreover, X and Y contain constraints that define the structure of the energy network. In particular, these constraints enforce energy productions, P t,s,l,p,c , that do not exist in the energy network to be zero. Since these constraints do not affect our analysis, we omit to report them and refer the reader to Babonneau et al. (2017) for a complete list of constraints.

A Demand Response Robust Version of ETEM
The ETEM model presented in Section 3 optimizes the future evolution of the energy system by assuming full information. However, in reality, there is always uncertainty. In particular, demand response is an important source of uncertainty. ETEM assumes that demand response is completely controlled by the supply side. From a decentralized electricity market point of view, this structure can be interpreted as a situation where the supply side targets an optimal level of DR, and encourages the actual demand load to be close to target DR using DR programs such as incentive real-time pricing. But, one must expect some realized deviations from the plans. Such deviations can either be caused by differences in the actual overall daily demand loads, as time evolves over a 50-year horizon, or by the level of participation and compliance of consumers in DR programs. In this paper, we model both types of uncertainties as implementation errors of the planned demand response V and we immunize ETEM against such error (see Ben-Tal et al. 2015, for a seminal work on robust optimization with implementation error). Specifically, given a set C U ⊆ C D of useful demand, such as residential electricity, with error-prone demand response, the demand Θ t,l,c V t,s,c for time period t ∈ T, time-slice s ∈ S, region l ∈ L, and commodity c ∈ C U , is replaced with Θ t,l,c V t,s,c + δ t,s,c , where δ t,s,c captures the Relative Demand Response Deviations (RDRD) from the planned response. In a general presentation, we will consider the vector of relative deviations δ ∈ R d , to be the vector that is obtained from arranging the parameters δ t,s,c for all t ∈ T, s ∈ S and c ∈ C U into a vector. For simplicity of exposure, we assume that the vector of each seasonal subset of RDRD lies in a scaled budgeted uncertainty set (as introduced in Bertsimas and Sim 2004) defined as: where d := |T| · |S| · |C U |, β t,s,c captures the maximum possible RDRD , while Γ t,j is a seasonal budget that represents the fact that we expect a maximum of Γ t,j RDRD to take on their extreme values in season j of time period t. The idea of decomposing the uncertainty information structure over each season has three important advantages. First, from a numerical perspective, it will help by enabling us to employ decomposition schemes based on Benders decomposition. Second, it is also attractive from a statistical perspective as it allows us to calibrate the size of each seasonal uncertainty set using seasonal data which is usually readily available, while data over joint realizations over many seasons can be scarce. Finally, imposing uncertainty budgets on subsets of perturbed parameters is known to lead to less conservative solutions in robust optimization as it enforces the worst-case scenario to distribute the damage over parameters of different subsets instead of focusing on a single In what follows, we present first, in Section 4.1, how we modify ETEM in order to identify capacity and demand response plans that are immunized against RDRD. The modified model is cast as a robust multi-period linear program. Unfortunately, this class of problems is known to be generally intractable. To address this issue, Section 4.2 presents a tractable conservative approximation that can be reformulated as a large-scale linear program.

A Robust Multi-Period Linear Programming Formulation
With the introduction of the RDRD in ETEM, the problem has the potential of turning into a multi-stage decision-making problem where different decisions are allowed to depend on the observed realization of RDRD . Indeed, we will assume that in the first stage the energy network planner decides on the capacity expansion of the system (variable C) and the planned demand response (variable V ). Then, in later stages, on a season-byseason basis, the seasonal RDRD are revealed to the planner who responds by deciding on the optimal energy production, energy transfers between regions, and energy imports and exports. To simplify the analysis of such a multi-stage model, we will further assume that the planner's response only depends on the current seasonal RDRD instead of reacting to all the observations made since the first season.
To be mathematically precise and help with presentation, decision variable x ∈ R m , with m = |T| · |L| · |P| + |T| · |S| · |C D |, will capture the vector of all first stage decisions (C, V ) in the demand response robust version of ETEM. Furthermore, for each time period t and season j, let i ∈ I denote the seasonal period (t, j), and consider that y i ∈ R ki , with k i = |S j | · |L| · |P| · |C| + |S j | · |L| · |C I | + |S j | · |L| · |C EX | + |S j | · |L| 2 · |C T R |, captures all adaptive "procurement" decisions (P t,s,l,p,c , I t,s,l,c∈C I , E t,s,l,c∈C EX , T t,s,l,l ,c∈C T R ) t,s∈Sj ,l∈L,l ∈L implemented in season j and time period t. For simplicity, we will refer to seasonal periods as i ∈ I instead of (t, j) ∈ T × J and consider A i := {(t, s)} s∈Sj as the set of all (t, s) pairs associated to seasonal period i. Finally, we refer to ζ i ∈ [−1, 1] di as the vector of seasonal RDRD perturbances {ζ t,s,c } (t,s)∈Ai,c∈C U .
The chronology of our multi-period decision-making model is depicted in Figure 3. In particular, this model can be shown to take the following form: where constraint (3c) describes the polyhedral feasible set for C and V based on X and constraints (1f) and (1g), while constraint (3b) describes for each i the constraints imposed on {P t,s,l,p,c , I t,s,l,c∈C I , E t,s,l,c∈C EX , T t,s,l,l ,c∈C T R } (t,s)∈Ai,(l,l )∈L 2 . Namely, the latter includes the constraints in Y, which in Babonneau et al. (2017) decompose over each season, and constraints (1b)-(1e), (1h)-(1m). Finally, models the linear effect of the standardized perturbances on each constraint. For completeness, we note that f ∈ R m , h i ∈ R ki , A i ∈ R ni×m , B i ∈ R ni×ki , and b i ∈ R ni .

Affinely Adjustable Approximation
Given the known numerical intractability of problem (3) (see Ben-Tal et al. 2004), one can instead solve a conservative approximation model obtained by enforcing that adjustable variables y i (ζ i ) be affine functions of ζ i ; i.e., y i (ζ i ) = y i + Y iζi , where y i ∈ R ki and Y i ∈ R ki×ri become the decision variables. We follow this approach after employing a commonly used lifting of the uncertainty set: and I di is the d i × d i identify matrix, r i = 2d i . In other words, problem (3) is conservatively approximated by the following robust linear program (referred as the Robust Multi-Period Conservative Approximation model): Applying standard robust reformulation techniques, we obtain the following tractable linear program: where ψ i ∈ R (di+1)×1 and Φ i ∈ R ni×(di+1) are new variables added to the model to guarantee the feasibility of the decision variables in the worst-case scenario. It is worth mentioning that problem (5) provides an upperbound and feasible first-stage decision for problem (3) as the space of decision rule has been reduced to affine functions.

Improving Numerical Efficiency using Benders Decomposition
Because of the capacity expansion decisions with a horizon of several decades, alongside with the procurement decisions with daily resolution, ETEM is already in its deterministic form a large-scale LP model. The size of LP that needs to be solved is further exacerbated when considering the robust reformulation presented in (5). In this regard, Section 6 will consider instances where the number of decisions and constraints is increased twenty-fold and nine-fold respectively. For this reason, we propose here a decomposition scheme to reduce the computational burden. Specifically, we start by reformulating problem (4) as follows: with where constraint (6d) is a redundant constraint that describes the fact that (4b) must be satisfied forζ i = 0 since 0 ∈Z i . Given that g i (x, y i ) is jointly convex, we identify a supporting plane representation that will enable the use of a constraint generation approach, also commonly referred to as a Benders decomposition scheme, which relies on iteratively solving a master and set of subproblems until convergence. Note that constraint (6c) should be considered violated if x and y i are such that problem (7) is infeasible.
In what follows, we start by describing the algorithm, then present two schemes that can be employed to improve its numerical efficiency by refining the quality of the upper and lower bounding process. In the latter case, the refinement will make use for the first time of the notion of Pareto robustly optimal solution (Iancu and Trichakis 2014) in decomposition schemes for robust optimization.
Remark 1. It is worth noting that our scheme is strongly inspired by the one used in Ardestani-Jaafari and Delage (2018) and Ardestani-Jaafari and Delage (2020), yet we depart from that work in three ways. First, unlike this prior work we choose to keep in problem (6) the y i variables, which describe what the affine strategy prescribes for the nominal scenarioζ = 0. This allows us to more easily interpret the optimal solution. We also observed empirically that this led to a significant reduction in the number of iterations. Second, we derive conditions that are specific to this application under which we can ensure that the solution of the master's problem always leads to a set of feasible subproblems. Third, we identify a new scheme, which is based on the theory of Pareto robust optimality, to identify optimal solutions of the master problem that accelerate the convergence of the algorithm.

The Benders Decomposition Algorithm
To simplify the presentation of this algorithm we make the following assumption. Assumption 1. (Uncapacitated technology) ETEM includes a technology p x ∈ P P c /P R with c ∈ C U generated from a resource c x ∈ C I /C D , with infinite capacity, i.e., Ω t,l,px = ∞. Furthermore, X and Y do not impose additional constraints on P t,s,l,px,c and I t,s,l,cx .
The role of the above assumption is to ensure that when removing constraint (6c) from problem (6), we do not end up with solutions for (x, {y i } |I| i=1 ) that make problem (7) infeasible. Note also that in our description of ETEM, the EXT technology and ExR energy resource respectively play the role of p x and c x .
Proof. One needs to show that for any given x and y i , i = 1, · · · , |I|, that satisfy (6b) and (6d), there exists a Y i so that the constraint (4b) is feasible. To do so, given that we already know x and y i are such that they can cover the nominal demand response, i.e., ζ + = ζ − = 0, we will show how to construct a Y i that will simply cover any deviation of the demand coming from the RDRD using the uncapacited technology. Specifically, we let all values of Y equal to zero except for the terms that model the influence of δ + on P t,s,l,px,c with c ∈ C U , P t,s,l,px,cx , and I t,s,l,cx . In particular, when setting all other terms to zero, most constraints reduce exactly to the constraints imposed in (6d) and are straightforwardly satisfied. The only constraints that need to be verified consist of a few members of (1b), (1e), (1h) and (1j), which can be described as follows: whereP ,V ,C,Ī,T , andĒ refer to the assignments made through the fixed x and y i 's, while P + t,s,l,px,c , P + t,s,l,px,cx,c , and I + t,s,l,cx,c model the influence of δ + t,s,c on P t,s,l,px,c , P t,s,l,px,cx , and I t,s,l,cc respectively, for c ∈ C U . In order to get a feasible assignment for Y , for all t, s, c ∈ C U , we let P + t,s,l,px,c := Θ t,l,c , P + t,s,l,px,cx,c := P + t,s,l,px,c /η t f,f , and I + t,s,l,cx,c = P + t,s,l,px,cx,c /η cx . One can readily verify that constraints (8), (10), and (11) reduce to the constraint already accounted for in (6d). On the other hand, constraint (9) reduces to: which is necessarily satisfied since Θ t,l,c ≥ 0 and δ − i ≥ 0. A similar argument can be used to confirm that all non-negative constraints onP t,s,l,px,c + P + t,s,l,px,c δ + t,s,c ,P t,s,l,px,cx + c∈C U P + t,s,l,px,cx,c δ + t,s,c , andĪ t,s,l,cx + c∈C U I + t,s,l,cx,c δ + t,s,c are also satisfied.
In practice, if Assumption 1 is not satisfied, a few remedies exist. First, one can create an artificial technology with either infinite or very large capacity, whose price is set to an arbitrarily large amount. One can interpret the price of this artificial technology as the marginal cost of shortages in the system. If this technology is not used in the optimal solution that is identified, then shortages do not occur and the solution is optimal for the problem without this artificial technology. Alternatively, one could also modify the algorithm that is presented below to have it identify "feasibility cuts" when the current candidate solution makes (7) infeasible (see Rahmaniani et al. 2017, for details).
We now focus on presenting a support plane representation of g(x, y i ).
Theorem 2. Given that Assumption 1 is satisfied, Proof. In the first step, we identify the robust counterpart of constraint (4b). For each i ∈ I, given thatZ i is non-empty, LP duality applies on each robust constraint thus introducing auxiliary variables that we denote in matrix form in Φ i . Namely, constraint (4b) can be said equivalent to Therefore, problem (7) can be reformulated as: which is a feasible minimization problem according to Lemma 1 since we assumed that Assumption 1 holds. In the next step, sinceZ i is a compact and convex set, Sion's minmax theorem (Sion 1958) holds and one can reverse the order of minimization over {Y i , Φ i } and maximization overζ i 's. Finally, since the set defined by (13) and (14) is feasible, one can employ LP duality to replace the inner minimization problem with its dual maximization problem to obtain: g(x, y i ) = max ζ i ∈Zi,θi,λi where λ i ∈ R ni×ri , θ i ∈ R ni×1 are dual variables associated with constraints (13) and (14) respectively.
Equipped with Theorem 2, we can define the following Benders decomposition algorithm. Intuitively, the procedure consists in solving a so-called master problem in which constraint (6c) is removed thus producing a sequence of lower bounds for problem (4). Given a current optimal solution of the master problem, the latter is refined by progressively reintroducing the constraints that are the most violated from the set: For each i, the most violated constraint can be found by solving the LP presented in (12), which optimal value can be used to update an upper bound problem (4). We refer the reader to the pseudo-code presented in Algorithm 1.
Algorithm 1: Benders Decomposition (BD) algorithm 2 Let i be the index for each subproblem; 3 Let v = v + 1 and solve the deterministic problem, i.e., problem (4) whenZ i = {0}, and store the value of x as x (v) and y i as y i ) and store the optimal value in ρ * i and decision variables in Solve the updated master problem : Let v = v + 1 and store the optimal value of the master problem variables in x v , y v i and ρ v ; 8 Update the lower bound:

Improving the Algorithm's Convergence
In order to improve the algorithm convergence, we apply two strategies. First, as proposed in Ardestani-Jaafari and Delage (2018) and Ardestani-Jaafari and Delage (2020), we add a set of valid inequalities to the master problem (15) based on violated scenarios from previous iterations in order to tighten the lower bounding problem. Specifically, for any i ∈ I and finite set of scenarios {ζ l i } l∈Ω ⊂Z i , the following constraint necessarily holds in problem (6): given that for all l ∈ Ω: where y l i plays the role of an adjustable plan specifically designed for scenarioζ l i and which can depart from the affine decision rule. Finally, In our implementation, we choose to set {ζ l i } l∈Ω to be the violated scenarios in the last v iterations of the algorithm. For completeness, we include the new master problem below: The second proposed improvement has to do with how the upper bound is obtained. In particular, this bound comes from evaluating the true worst-case cost of the current best candidate based on the master problem (15), or its tighter version (16). Yet, since the master problem is itself a robust optimization problem, we can expect based on the findings of Iancu and Trichakis (2014) that at each iteration there exists a large set of optimal solutions for (15). In a traditional Benders decomposition approach (such as in Ardestani-Jaafari and Delage 2018), the candidate that is used to calculate the upper bound and generate an optimality cut is arbitrarily chosen by the LP solver used to solve (15). While such a solution is optimal for (15), it could have a much worst performance in problem (12) compared to other optimal solutions. We therefore recommend, after solving (15), to identify a Pareto robustly optimal solution of (15) by solving the following LP: where M * v is the optimal solution of the master problem (15) in iteration v of Algorithm 1, and > 0 allows for some small sub-optimality in order to help numerically. The new candidate (x v , y v i ) remains approximately optimal for (15) yet is guaranteed to not be Pareto dominated by other solutions of (15). This follows from the fact that constraint (17c) is equivalent to: and that ( Corollary 1 in Iancu and Trichakis (2014) can be applied to obtain the Pareto non-dominance guarantee. We refer the reader to Iancu and Trichakis (2014) for more details about how to identify Pareto robustly optimal solutions.

Computational Results
In this section, we evaluate the performance of the proposed robust ETEM approximation model (i.e., RMPCA) on a case study based on the energy system of the Arc Lémanique region in Switzerland. To do so, we compare the performance of the RMPCA formulation with a static robust model (SRM), which considers all decisions as here-and-now decisions, and a deterministic model (DET), which disregards demand response uncertainty. The purpose of these computational experiments is to empirically show that: i) considering DRU in a capacity expansion planning problem and solving it using robust optimization can decrease both worst-case and expected total costs of the system compared to a deterministic formulation of the problem; and ii) assuming flexibility of the procurement decisions with respect to observed actual DR decreases the level of conservatism of the model and consequently the expected total cost of the system. Recall that in RMPCA the planner first decides on the capacities and planned DR, then in seasonal periods, after observing the actual DR, he decides on the optimal procurement of energy. Theoretically, it is known that the policies proposed by SRM are more conservative compared to policies proposed by RMPCA. In addition, we present a detailed comparison of the structure of the RMPCA, SRM, and DET policies. The purpose is to understand how considering the uncertainty of the demand response affects the long-term capacity expansion strategies. Finally, we provide a computational analysis to compare the performance of the different decomposition algorithms presented in Section 5. All numerical studies are performed on a 64-bit computer with 128 GB of RAM and Intel(R) Xeon(R) 3.1 GHz (40 CPUs). The deterministic version of ETEM is originally written in AMPL. However, we extract the standard form matrices from AMPL and carry out robust optimization using the YALMIP toolbox and GUROBI solver in a MATLAB R2019a environment.

A Swiss Case Study
The ETEM model describing the Arc Lémanique region (Cantons of Geneva and Vaud, in Switzerland) encompasses 142 different technologies, including centralized electricity and heat production plants, decentralized electricity production, conventional and flexible loads, and end-use transportation technologies. In addition, 57 energy commodities, including 21 types of useful demand, are modeled. We choose to immunize ETEM against residential electricity demand response deviations (i.e., Res-ELC in Figure 1), which implies that the robust model is subject to 120 independent sources of perturbations. This choice is motivated by the fact that the response behavior of residential electricity is more prone to uncertainty in the long-run compared to industrial and large-scale electricity consumers, as the former corresponds to a larger and more diverse group of consumers. Table 1 presents techno-economic features of both presently available and expected future electricity generation technologies in the Arc Lémanique region. As we wish to avoid as much as possible energy shortages in the system, we set the cost of the additional energy source (ExR) high enough so that there is no capacity shortage in the optimal solution of RMPCA model when the largest amount of deviations are assumed. Finally, in order to better investigate the influence of robust optimization on the solutions, no upper bound limits the installation of renewable technologies (wind and photovoltaic power plants) in the model. For interested readers, we note that our deterministic version of ETEM presents 171,620 variables and 431,652 constraints while the LP reformulation of RMPCA, i.e., problem (5), presents 3,540,720 variables and 3,708,658 constraints. While we refer the reader to Babonneau et al. (2017) for a complete description of all input parameters and assumptions, it is worth mentioning that we consider here a case where CO 2 emissions are curbed to 2.5 Mt in 2050, a 45% reduction compared to 2010 levels (5.48 Mt).

Performance Analysis
In this section, we study the performance of strategies proposed by RMPCA, SRM, and DET in the average and worst-case scenarios. For a specific budget size Γ and box size β, x Γ,β RM P CA represents the RMPCA strategy obtained by solving problem (5). Similarly, x Γ,β SRM is the static robust strategy obtained by solving problem (5) when forcing Y i = 0 ∀i ∈ {1 · · · |I|}. It should be noted that, because problem (5) is computationally demanding and cannot be solved in a 72-hour time limit, we have used the PRO-BD-VI algorithm proposed in Section 5. Finally, the DET strategy x DET is obtained by solving the deterministic model, i.e., problem (4) whenZ i = {0} ∀i ∈ {1 · · · |I|}. After obtaining all the strategies from the different models, in a next step, we generate a set of 1000 DR random scenarios ζ j i (∀i ∈ {1, · · · , |I|}, j ∈ [1, · · · , 1000]) using a uniform distribution on the [−β, β] interval. For each scenario, we then re-optimize problem (4) with the value of x and ζ i fixed to the corresponding strategy and random scenario. In the first period ( i.e., t = 1), we have no uncertainty given that a full information about the level of DR is available. If, for a specific policy and scenario, problem (4) becomes infeasible, this implies a capacity shortage. Table 2 reports the proportion of scenarios with capacity shortage for the different strategies and different budget sizes (expressed in percentage of |d i | and denoted as Γ%). When the budget size increases, the chances of capacity shortage of RMPCA and SRM reduces. This is because a larger budget size increases the conservatism of the model, and consequently, the model installs more capacity to avoid shortage, which obviously is more expensive. Nevertheless, for budget sizes in the range of [50%, 60%], the capacity shortage of RMPCA can be as low as 2%. To have a fair comparison between the expected total cost of different strategies, we report the average total cost both over all scenarios and over the scenarios without capacity shortage.  Figure 4a presents statistics of the total cost for the RMPCA and SRM strategies, over all scenarios, for different budget sizes and β = 0.6. Figure 4b gives similar results when the average costs are taken over only scenarios without capacity shortage. We can first note that the average cost for the SRM strategy remains fixed for budget sizes greater than 25%. This is because when Γ% = 25%, |d i | = 1, the SRM already considers each constraint to be maximally perturbed due to the fact that each one only involves one demand deviation. Second, we note that Figures 4 (a) and (b) show that the average RMPCA cost is considerably lower than the SRM cost for all budget sizes. This already suggests that the flexibility of energy procurement decisions can significantly reduce the expected total cost. More precisely, the maximum difference occurs when Γ% = 30% yet this difference comes at the price of 29% chances of capacity shortage (see again Table 2). A better trade-off is achieved when Γ% = 50%, where the average total cost of RMPCA is around 4% less than for SRM (i.e., a 9 billion CHF reduction) while the chances of capacity shortage are estimated at only 4% . In comparison, the average total cost (including shortage penalties) from implementing the DET strategy is around 256 billion CHF, which is about 33% more than the total cost of RMPCA. Overall, this evidence supports the importance of: i) accounting of uncertainty of the demand response in ETEM; and ii) the importance of considering that the procurement is adjustable with respect to the actual DR. Figure 5 is similar to Figure 4, but we have fixed the budget size to 50% and let the box size β vary. For example, β = 0.2 means that the actual demand response at each time-slice can deviate from the planned demand response by up to 20% of the total demand. As expected, increasing the support of ζ leads to an increase of the average total cost of the system. Interestingly, the rate of increase of the DET cost is larger than for RMPCA and SRM. Indeed, as SRM and RMPCA are sensitive to worst-case values of DR, more capacity is built under these strategies. This gives more production flexibility with cheaper prices for the system when the actual DR is lower than planned DR. Finally, the figure also illustrates that the average total cost (∼192.5 billion CHF) of preventing 96% of shortages when β = 0.6 is only about 5% more expensive than the cost estimated by the deterministic model (i.e. ∼182.5 billion CHF). Figure 6 corresponds to the Cumulative Distribution Function (CDF) diagram of the simulated total cost of RMPCA, SRM and DET strategies when Γ% = 50% and β = 0.6. We can note that RMPCA stochastically dominates SRM, while SRM stochastically dominates DET. In fact, the maximum cost of the RMPCA strategy over all random scenarios is less than the minimum cost of the SRM strategy, and similarly for the DET strategy.  Figure 4: Statistics of the total cost of RMPCA and SRM strategies for different budget sizes and β = 0.6. (a) presents the distribution of the total cost over all 1000 random scenarios, whereas (b) presents the conditional distributions given that no shortages occurred. The curves present the mean values while each box presents the min, max, 25th and 75th percentile and some extreme realizations (+).   Again, this figure shows the importance of considering uncertainty of the DR in capacity expansion models. Next, Figure 7 shows the worst-case total cost of the system when implementing the RMPCA and SRM strategies for different budget sizes. As this figure shows, RMPCA performs up to 6% better than SRM which is equivalent to approximately 12 billion CHF over the entire planning horizon. Figure 8 provides finally an analysis on the trade-off between rate of capacity shortage and the average total cost for RMPCA and SRM strategies when the budget size varies from 10% to 100% and the box size β = 0.6. As it is also shown in Table 2, for budget sizes greater than 30%, the rate of capacity shortage for SRM is 0%. However, this zero-shortage comes at the cost of building too much capacity (see Figure 9, below). As Figure 8 suggests, allowing a certain degree of shortage can reduce the total cost drastically. For example, for Γ% = 60% allowing a 2% chance of capacity shortage reduces the average total cost by around 7 billion CHF (equivalent to 4% of the total system cost).

Electricity Generation Structure
In this section, we discuss how the different modeling schemes impact the evolution over time of electricity generation. Figure 9 presents first the total installed capacity for the different strategies. SRM is the most conservative approach, as it yields more capacity than RMPCA and DET (around 1.5 times more than RMPCA, and 3 times more than DET). It is indeed reasonable that RMPCA is less conservative than SRM, as it enjoys the flexibility of the energy procurement decisions to adjust for the actual (observed) demand response. Figure 10 presents next the installed capacity by type of technology for the RMPCA, SRM and DET strategies. It highlights the importance of wind power to decarbonize the electricity generation sector. Specifically, wind power will form around 88% of the total installed capacity by 2050. It is followed by hydro power plant with 8%, and other technologies (including CHP, PV, geothermal, fuel-cell, and municipal-waste power plants) with less than 4% share in the total installed capacity. In addition, Figure 11 reports on the average annual electricity generation by different technology types for RMPCA, SRM, and DET policies. A more detailed description of the installed capacity and the average annual electricity generation is presented in Appendix B.

Computational Analysis
In this section, we numerically evaluate the performance of the decomposition algorithms presented in Section 5. Specifically, we investigate how considering PRO and valid inequalities can improve the solution time and the number of iterations of Algorithm 1. In order for our results to be more general, we have designed smaller versions of ETEM with 16 technologies and 9 energy commodities, in which parameters are randomly assigned. Moreover, the size of the model changes with different problem horizons T ∈ {3, 5, 7, 9, 11, 13}. Then, for each problem size, we generate 20 random instances and solve these random instances for different budget sizes Γ% ∈ {20%, 40%, 60%, 80%}. Table 3 reports, for different problem sizes, the average solution time (ST) of CPLEX for solving the RMPCA problem (5), together with the average solution time and the number of iterations (It) for different versions of our decomposition algorithm. Specifically, BD refers to the traditional Benders decomposition, i.e., Algorithm 1. PRO-BD refers to the algorithm with improved upper bound based on PRO solutions while PRO-BD-VI additionally exploits valid inequalities. It is worth mentioning that based on our numerical experiments, restricting Ω to the last 4 worst-case scenarios led to the best solution time.
Results show that PRO-BD improves the solution time of the BD by up to 36%, whereas PRO-BD-VI is able to improve it by up to 59%. On average, PRO-BD improves the solution time by 22% and PRO-BD-VI improves by 40%. Finally, a considerable improvement in terms of the number of iterations is observed when comparing PRO-BD-VI and BD.

Concluding Remarks
In this paper, we model demand response uncertainty in an adjustable multi-period robust generation capacity expansion problem. In a first stage, the planner decides on capacity expansion, as well as, planned demand response. Then, in sequence of seasonal stages, the actual DR is revealed, and accordingly, the optimal procurement decisions are taken. As the resulting adjustable model is intractable, we use a conservative approximation scheme based on affine decision rules to solve the problem. We present a general formulation of the problem, where with slight modifications of the formulation one could accommodate other sources of uncertainty such as the level of production one can achieve from intermittent resources. We test our algorithm on a real-world case study based on the energy system of the Arc Lémanique region (Cantons of Geneva and Vaud, in Switzerland). Numerical results reinforce the importance of adjustable robust formulation. Namely, RMPCA can reduce the average total cost of the system by 33% compared to the solution of the deterministic problem, while keeping the chances of shortage near 4%. RMPCA also performs 4% better than a more naive static robust counterpart (SRM). Finally, we develop a Benders type of decomposition to solve our large-scale RMPCA problem. In this algorithm, the speed of convergence is improved through i) adding valid inequalities to the master problem, and ii) identifying PRO solutions of the master problem. One alternative strategy to calculate the optimal investment decisions is to sequentially solve forward-shifting horizon RMPCA problems, and provide the investment plan over the entire horizon by assembling the capacity addition decisions made at the beginning of each forward-shifting horizon. However, this practical implementation of RMPCA problem is beyond the scope of this paper and shall be considered as future work. Moreover, it is worth mentioning that although this paper has focused on the uncertainty of demand response, one can easily introduce a group of supply-side uncertainties, such as availability of intermittent resources, or uncertainty of investment costs, to the general formulation of the problem (4). In this case, matrices A i and vector f will be affected by these uncertainties. As our problem still preserves the "fixed recourse" condition, the approach presented in Sections 4.2 and 5 remains valid. Conversely, considering another group of supply-side uncertainties, such as fuel cost, or efficiency uncertainties, will affect matrices h i and B i , and thus consists of more significant changes to the structure of the decision model. In conclusion, here follows some interesting directions to expand on this paper. One idea is to consider investment decisions to be adaptive to the evolution of the uncertainty. Another extension is to consider other sources of uncertainties, especially the ones that do not preserve the fixed recourse characteristic of the problem. One could also control the loss of loads through reliability constraints. Finally, another possible extension would be to link the robust ETEM with a collective behavior demand model to more precisely adjust the size of the specific uncertainty sets of each period.