Samuel Torres-Rincón

Mauricio Sánchez-Silva

Emilio Bastidas-Arteaga

A Multistage Stochastic Program for the Design and Management of Flexible Infrastructure Networks

Keywords: Flexibility, Multistage Stochastic Programming, Value of Flexibility, Infrastructure Networks

Throughout their lifetime, infrastructure network systems face unplanned events that impose pressures on their integrity, functionality, and ability to deliver value. Most of the existing infrastructure is designed to deal with the challenges imposed by uncertain external phenomena. Authors from different backgrounds have identified flexibility, changeability, and adaptability as key attributes that modern systems should have to face uncertain scenarios. Specifically, flexibility is an ability that allows a system to be easily adapted when necessary. The concept of flexibility is compelling, but it is not clear how to measure the value it may provide. Determining how much to pay to introduce flexibility is an essential aspect of designing flexible systems, but the dependence of this value on the future evolution of the system results in a complex decision process. The sequential nature of the process can be modeled using multistage stochastic programming. The model explicitly considers the flexibility built into the network components as a decision variable at the initial stage. The model is tested in a generic infrastructure network that must meet a stochastic demand.

The results show the relationship between the value of flexibility and the life-cycle costs at the construction and operation stages.

Introduction

Background

Infrastructure networks, as any infrastructure system, are constantly exposed to ever-changing conditions. Changes in the demand, the supply sources (in the case of water and energy distribution systems), and regulations may impose large stresses on the system [START_REF] Browning | Measuring the Life-Cycle Value of Enduring Systems[END_REF][START_REF] Vieira | Systemic approach for the capacity expansion of multisource water-supply systems under uncertainty[END_REF][START_REF] Zhang | Resilience-based network design under uncertainty[END_REF][START_REF] Sánchez-Silva | Managing Infrastructure Systems through Changeability[END_REF]. This uncertainty can either challenge or boost the ability of the network to perform at required levels [START_REF] Lozano | Developing flexible management strategies in infrastructure: The sequential expansion problem for infrastructure analysis (SEPIA)[END_REF]. In the former case, networks usually rely on the robustness provided by over-dimensioned designs which can be largely inefficient, as in the case of sanitation systems [START_REF] Roefs | Centralised, decentralised or hybrid sanitation systems? Economic evaluation under urban development uncertainty and phased expansion[END_REF]. In the latter case, even if the new conditions are favorable, the network may not have the tools to take advantage of the new opportunities [START_REF] Cardin | Improving the Lifecycle Performance of Engineering Projects with Flexible Strategies: Example of On-Shore LNG Production Design[END_REF]. To face these challenges, many authors have identified the need to develop flexible and adaptable systems, which have the potential to improve the sustainability and efficiency of the system under highly uncertain conditions [START_REF] De Neufville | Flexibility in Engineering Design[END_REF][START_REF] Mortazavi-Naeini | Application of multiobjective optimization to scheduling capacity expansion of urban water resource systems[END_REF][START_REF] Fraga | Planning for infrastructure capacity expansion of urban water supply portfolios with an integrated simulation-optimization approach[END_REF][START_REF] Erfani | Real-Options Water Supply Planning: Multistage Scenario Trees for Adaptive and Flexible Capacity Expansion Under Probabilistic Climate Change Uncertainty[END_REF][START_REF] Chester | Toward adaptive infrastructure: flexibility and agility in a non-stationarity age[END_REF]. Specifically, flexible systems are defined as systems with the ability to change as "easily" as possible [START_REF] Ross | Defining Changeability : Reconciling Flexibility , Adaptability , Scalability , Modifiability , and Robustness for Maintaining System Lifecycle Value[END_REF][START_REF] Saleh | Flexibility : a multi-disciplinary literature review and a research[END_REF][START_REF] Špačková | Long-term adaption decisions via fully and partially observable Markov decision processes[END_REF][START_REF] Sánchez-Silva | Managing Infrastructure Systems through Changeability[END_REF][START_REF] Torres-Rincón | Conceptual and Numerical Analysis of Flexibility in Infrastructure Systems[END_REF]. This can be measured as the money or time investment required to modify the system. In the case of infrastructure networks, adaptations can be usually seen as the addition (or removal) of nodes and links or the expansion (or downgrade) of their capacity. Other changes, mostly managerial, may be possible (e.g. reversal of the flow direction). Flexibility is introduced by investing a number of resources during the design and construction phase to enable the option of deploying future adaptations at lower costs and shorter lead times. This investment is usually called the value of flexibility and is defined as the amount the system stakeholders' are willing to pay to introduce flexibility into the system [START_REF] Cardin | Improving the Lifecycle Performance of Engineering Projects with Flexible Strategies: Example of On-Shore LNG Production Design[END_REF][START_REF] De Neufville | Flexibility in Engineering Design[END_REF]. It is assumed that the enabled adaptations are capped at a certain value that depends on the amount invested at the initial stages of the project [START_REF] Torres-Rincón | Conceptual and Numerical Analysis of Flexibility in Infrastructure Systems[END_REF]. This interpretation of flexibility is similar to the concept of shell capacity described by Angelus et al. [START_REF] Angelus | Optimal Sizing and Timing of Modular Capacity Expansions[END_REF].

Previous works and applications

While flexibility is not limited to capacity expansion, this type of problem is the most common and has been extensively studied in the manufacturing sector [START_REF] Luss | Operations Research and Capacity Expansion Problems: A Survey[END_REF][START_REF] Van Mieghem | Capacity Management, Investment, and Hedging: Review and Recent Developments[END_REF][START_REF] Martínez-Costa | A review of mathematical programming models for strategic capacity planning in manufacturing[END_REF]. The first work comes from Manne [START_REF] Manne | Capacity Expansion and Probabilistic Growth[END_REF], who modeled the optimal excess capacity of pipelines, highways, and steel plants subject to a random walk demand. Luss [START_REF] Luss | Operations Research and Capacity Expansion Problems: A Survey[END_REF] further detailed the problem as determining the size, location, and timing of future adaptations while minimizing the discounted costs of the expansion processes. This author also identified the modeling of the demand, the representation of the capacity as a continuous or discrete variable, the selection of the discount rate, and the presence of economies of scale as key elements of the problem formulation. Mathematical programming rapidly became one of the preferred methods to solve this type of problem, as shown by the work of Rajagopalan [START_REF] Rajagopalan | Capacity expansion with alternative technology choices[END_REF], who formulated a mathematical programming model to determine the optimal initial capacity and the technology acquisition decisions, for increasing demands. Later came the work of Ahmed and Sahinidis [START_REF] Ahmed | An Approximation Scheme for Stochastic Integer Programs Arising in Capacity Expansion[END_REF] who developed a multistage stochastic (MS) integer program to analyze a multiperiod investment model for capacity expansion under uncertain demand and costs, and including the effect of economies of scale. Singh et al. [START_REF] Singh | Dantzig-Wolfe Decomposition for Solving Multistage Stochastic Capacity-Planning Problems[END_REF] also developed an MS integer programming model for planning discrete capacity expansion of production facilities. Similarly, Huang and Ahmed [START_REF] Huang | The Value of Multistage Stochastic Programming in Capacity Planning Under Uncertainty[END_REF] formulated a general MS capacity planning model with discrete capacity for the semiconductor industry.

The prevalence of stochastic and dynamic programming models to solve capacity expansion problems in the manufacturing sector led to the adoption of these techniques in other areas such as planning and design of infrastructure networks. For instance, in the context of transportation networks, Marín and Jaramillo [START_REF] Marín | Urban rapid transit network capacity expansion[END_REF] formulated a multi-period capacity expansion problem for rapid transit network design. Karoonsoontawong and Waller [START_REF] Karoonsoontawong | Integrated Network Capacity Expansion and Traffic Signal Optimization Problem: Robust Bi-level Dynamic Formulation[END_REF] developed a robust optimization model for the problem of continuous traffic network capacity expansion with dynamic traffic assignment and traffic signal optimization. Gao et al. [START_REF] Gao | Integrated maintenance and expansion planning for transportation network infrastructure[END_REF] considered the case where a road network is expanded either by adding new links or by increasing the capacity of the existing links, integrated with the road maintenance problem using a mixed-integer, non-linear, bi-level optimization program, using multi-period decisions. In the area of water and sewage distribution systems, Mortazavi-Naeini et al. [START_REF] Mortazavi-Naeini | Application of multiobjective optimization to scheduling capacity expansion of urban water resource systems[END_REF] developed a multi-objective optimization approach for the planning of urban water system expansions that considers the combined effect of operating rules and infrastructure conditions. Saif and Almansoori [START_REF] Saif | A capacity expansion planning model for integrated water desalination and power supply chain problem[END_REF] developed a model for the expansion of water desalination and power supply infrastructure using a deterministic multi-period mixed-integer linear programming formulation.

Similarly, Fraga et al. [START_REF] Fraga | Planning for infrastructure capacity expansion of urban water supply portfolios with an integrated simulation-optimization approach[END_REF] used dynamic programming to develop an integrated framework for the optimization of water supply system expansions, considering short and long-term water supply sources. In the case of energy production and transmission networks, Loureiro et al. [START_REF] Loureiro | Capacity expansion in transmission networks using portfolios of real options[END_REF] used the concept of real options combined with a mixed-integer linear programming model for multistage expansion planning. Cardin et al. [START_REF] Cardin | Strategic real option and flexibility analysis for nuclear power plants considering uncertainty in electricity demand and public acceptance[END_REF] combined decision rules and stochastic programming to model flexibility in a nuclear infrastructure system, considering a random demand and including the social acceptance of nuclear power as a limitation to the system expansion capabilities. For networks in general, Taghavi and Huang [START_REF] Taghavi | A Multi-Stage Stochastic Programming Approach for Network Capacity Expansion with Multiple Sources of Capacity[END_REF] applied the concepts of spot market and contract capacity combined with an MS integer program to model a network with multiple sources of capacity.

In both manufacturing systems and infrastructure networks, the adaptation process requires large capital investments (in proportion with the total system size), which are sometimes irreversible [START_REF] Berman | A stochastic optimization model for planning capacity expansion in a service industry under uncertain demand[END_REF], and are susceptible to the negative consequences of large lead times. In both cases, there is a contradictory effect from economies of scale, which may favor both large initial designs and large future adaptations. However, infrastructure networks are subject to more strict governmental regulations, where failure to comply with minimum levels of performance can result in heavy societal costs. This poses an additional problem of strategically selecting the extent to which flexibility should be incorporated. Besides, infrastructure networks are usually planned for longer service horizons, which favors delayed deployments due to discounting [START_REF] Ryan | Capacity Expansion for Random Exponential Demand Growth with Lead Times[END_REF]. Furthermore, the nature of stakeholders affects the management and adaptation strategies. While manufacturing systems are generally completely privately owned, infrastructure network ownership can range from completely public to completely private, which may create conflicting interests. For instance, in energy infrastructure, maintaining a determined level of excess capacity may be valuable for society but not for profit maximization [START_REF] Garcia | Equilibrium capacity expansion under stochastic demand[END_REF]. Finally, the frequency of adaptation is considerably different: in manufacturing systems is not uncommon to see changes every one or two years, while infrastructure networks may see changes in periods of five to ten years.

Objectives and scope

The works discussed in the previous section show the potential of (and the preference for) MS programming as a tool to model and solve complex sequential decision problems. Considering that the problem of designing and managing flexible infrastructure networks is in itself a sequence of decisions under uncertainty, this paper proposes a multistage stochastic program (MSP) to analyze the problem of flexibility. Operational and maintenance costs are explicitly considered, to take into account the relative importance of these costs in the systems modeled and the impact these costs may have on the preference of the model for flexible solutions. The model uses scenario trees generated by Monte-Carlo simulation and k-medoids clustering to simulate and discretize the random process. The novelty of the model is that it explicitly considers the flexibility range introduced by design as a decision variable, which will restrict all the future planned adaptations.

The proposed model results in an alternative methodology to real options as a tool to determine the value of flexibility. By modeling flexibility as a decision variable, the program decides how much should be paid at the initial stage to lower future expenses.

The article is organized as follows: 2 gives an overview of the concept of flexibility in infrastructure systems and its numerical representation. A review of the theory behind two-stage and MSPs is presented in 3. 4 provides the formulation of the proposed MSP. 5 presents a methodology to solve the MSP using an approximate approach. Finally, ?? provide a numerical example with a generic network and a summary of the key results.

Designing for Flexibility

The concept of flexibility is considered to be not "academically mature" [START_REF] Saleh | Flexibility : a multi-disciplinary literature review and a research[END_REF] due to the lack of a precise and universally accepted definition. In this study, flexibility will be understood as defined by [START_REF] Cardin | Improving the Lifecycle Performance of Engineering Projects with Flexible Strategies: Example of On-Shore LNG Production Design[END_REF][START_REF] Špačková | Long-term adaption decisions via fully and partially observable Markov decision processes[END_REF], and [START_REF] Torres-Rincón | Conceptual and Numerical Analysis of Flexibility in Infrastructure Systems[END_REF]. The authors define flexibility as the ability of a system to easily adapt any of its components (or subsystems). The effort necessary to complete an adaptation is measured as the number of required resources, which are usually represented as a monetary quantity. Furthermore, any adaptation will be limited to a maximum value given by the system's context and the resources invested to make the system flexible. These additional costs may be incurred at the initial stage of the project to introduce the option to adapt some of the system components. These expenses may come from research and development activities or by installing physical elements that will facilitate future adaptations. For instance, the floating platforms used in offshore wind farms could be designed to be expanded or be built larger-than-required to facilitate future adaptations to the turbine sizes.

The problem of determining how much should be invested to have flexibility in the system is usually known as measuring the value of flexibility. This is a complex problem because capturing the value of having the option to modify the system under certain circumstances is not a straightforward process [START_REF] De Neufville | Flexibility in Engineering Design[END_REF]. Different approaches have been developed based on different assumptions and modeling frameworks. For instance, the Real Options Analysis (ROA) method was developed by adapting the concept of options in financial markets to physical assets [START_REF] Trigeorgis | Real Options and Interactions with Financial Flexibility[END_REF][START_REF] Deneufville | Real Options: Dealing With Uncertainty in Systems Planning and Design[END_REF]. Under this approach, the value of flexibility is measured as the pricing of an option that provides the right but not the obligation to modify the system in some predefined way. In a different approach, Cardin et al. [START_REF] Cardin | Improving the Lifecycle Performance of Engineering Projects with Flexible Strategies: Example of On-Shore LNG Production Design[END_REF] defined the value of flexibility as the difference between the expected net present value (ENPV) of a flexible system and an inflexible system subjected to the same conditions. Similarly, Špačková and Straub [START_REF] Špačková | Long-term adaption decisions via fully and partially observable Markov decision processes[END_REF] defined it as the additional investment that should be made to have a flexible system compared with the inflexible alternative, using Markov Decision Processes (MDP) to model the sequential decision process. Following these works, the value of flexibility will be understood in this paper as the maximum value that should be paid at the start of the system life-cycle to minimize the discounted costs of future adaptations and maximize the received utility.

In the area of chemical processes, the concept of Flexibility Analysis (FA) has been proposed to optimize the design and operation of chemical plants under uncertainty [START_REF] Swaney | An Index for Operational Flexibility in Chemical Process Design. Part 2: Computational Algorithms[END_REF][START_REF] Grossmann | Evolution of concepts and models for quantifying resiliency and flexibility of chemical processes[END_REF][START_REF] Ahmed | An Approximation Scheme for Stochastic Integer Programs Arising in Capacity Expansion[END_REF][START_REF] Zhang | On the relation between flexibility analysis and robust optimization for linear systems[END_REF]. The main objective of this approach is to determine the optimal design and control variables that allow feasible operation for the whole range of uncertain parameters. If feasible operation is achieved, then the process is considered flexible enough. To guarantee feasibility, FA uses the worst-case approach [START_REF] Zhang | On the relation between flexibility analysis and robust optimization for linear systems[END_REF] as in the case of robust optimization. In fact, FA and robust optimization share many concepts and methods, even if historically they evolved separately. One key difference, however, is the use of recourse variables (the control variables) to modify the plant's response to particular realizations of the uncertain parameters. While recourse variables are traditionally not used in robust optimization, they are a key concept in stochastic optimization. By combining concepts from these optimization approaches, FA can be used to determine the optimal design characteristics of flexible chemical plants.

The formulation of the flexibility concept and the optimization models proposed in FA may suggest that the approach can be extended to other engineering systems. Indeed, the design and management of flexible infrastructure networks can be properly modeled using the representation of design, state, and control variables from the FA approach. The recourse options provided via control variables can represent adaptation decisions instead. In both cases, flexibility exists to face the uncertainty in external phenomena. Despite these similarities, the concept of flexibility in infrastructure systems (and networks) is more closely related to the concept of Real Options [START_REF] Trigeorgis | Real Options and Interactions with Financial Flexibility[END_REF][START_REF] Deneufville | Real Options: Dealing With Uncertainty in Systems Planning and Design[END_REF], while the definition of flexibility in FA is more similar to the concept of robustness in infrastructure systems (as defined in [START_REF] Ross | Defining Changeability : Reconciling Flexibility , Adaptability , Scalability , Modifiability , and Robustness for Maintaining System Lifecycle Value[END_REF]). The problem of how much should be paid today to have the option to modify the system in the future is a very important notion in flexibility for infrastructure systems that does not exist in FA. Furthermore, the changes implemented using control variables in chemical processes are usually reversible, while modifying infrastructure systems involves a sequential decision process that results in incremental modifications. The sequential nature of the decision process requires a modeling approach that considers the temporal interdependence of the decisions (such as dynamic or stochastic programming) for which the tools developed in FA are insufficient. For these reasons and despite the conceptual similarities, the FA approach is not used in this paper and rather an MS programming model is proposed.

Numerical representation of flexibility

To formulate and solve the problem of the value of flexibility, it is useful to have a numerical representation of flexibility. Unfortunately, the lack of consensus in the conceptual definition is also present in the formal description. Various authors have proposed an array of indices with different levels of complexity [START_REF] Swaney | An Index for Operational Flexibility in Chemical Process Design. Part 2: Computational Algorithms[END_REF][START_REF] Ross | Defining Changeability : Reconciling Flexibility , Adaptability , Scalability , Modifiability , and Robustness for Maintaining System Lifecycle Value[END_REF][START_REF] Fitzgerald | Managing Uncertainty in Systems with a Valuation Approach for Strategic Changeability[END_REF][START_REF] Špačková | Risk based Optimization of Adaptable Protection Measures Against Natural Hazards[END_REF] but none has been universally adopted. This study will measure flexibility using the flexibility vector defined by Torres-Rincón et al. [START_REF] Torres-Rincón | Conceptual and Numerical Analysis of Flexibility in Infrastructure Systems[END_REF], shown in Equation 1. This formulation explicitly represents two characteristics usually associated with the concept flexibility: the presence of the option to change, and a measurement of the effort required to complete the change. The flexibility vector for the design or operation variable i has two dimensions: the first dimension measures the effort necessary to complete an adaptation as the ratio between the unitary cost of modifying the component without flexibility (without being specifically designed to be adapted) and the unitary cost of performing an adaptation when flexibility was introduced c nf,i /c f,i ; the second dimension measures the size of the available adaptation space as the ratio between the maximum value the design or operation property can take divided by its initial value x max,i /x 0,i . Therefore, a system can increase its flexibility by increasing the number of states that can reach through low-cost adaptations, or by reducing the cost of such adaptations.

fv i,t = c nf,i -c f,i c f,i , x max,t,i -x 0,i x 0,i (1)
For instance, Figure 1 shows two flexibility vectors of the same magnitude. Vector A, however, has a larger cost component (horizontal axis), while vector B has a larger adaptation space component (vertical axis). Clearly, the expected behavior of both designs under the same external conditions will be different. The vectorial representation of flexibility allows capturing the complex nature of the property.

Managing flexibility through policies

The elements of flexibility described previously are deeply interconnected with the physical characteristics of the system. A system is specifically designed and built to be flexible, which affects its morphology. Therefore, future management process to decide the optimal timing and magnitude of adaptations will be constrained by the limitations of the design [START_REF] Ross | Defining Changeability : Reconciling Flexibility , Adaptability , Scalability , Modifiability , and Robustness for Maintaining System Lifecycle Value[END_REF][START_REF] Fitzgerald | Managing Uncertainty in Systems with a Valuation Approach for Strategic Changeability[END_REF][START_REF] Cardin | Improving the Lifecycle Performance of Engineering Projects with Flexible Strategies: Example of On-Shore LNG Production Design[END_REF][START_REF] Cardin | Strategic real option and flexibility analysis for nuclear power plants considering uncertainty in electricity demand and public acceptance[END_REF][START_REF] Sánchez-Silva | Managing Infrastructure Systems through Changeability[END_REF][START_REF] Torres-Rincón | Conceptual and Numerical Analysis of Flexibility in Infrastructure Systems[END_REF].

The management process is also conditioned by external factors and the system's current state, but also on particular preferences of the stakeholders, technical limitations, user requirements, and regulatory frameworks. For instance, risk-averse stakeholders may prefer to build large systems from the beginning with enough flexibility to perform small adaptations, while stakeholders with a higher risk tolerance may favor smaller initial systems with high flexibility to deploy large adaptations.

These complex interactions between external elements and individual preferences can be modeled in the form of policies.

Policies are functions that map a set of states to decisions [START_REF] Cardin | Strategic real option and flexibility analysis for nuclear power plants considering uncertainty in electricity demand and public acceptance[END_REF][START_REF] Powell | Approximate Dynamic Programming: Solving the Curses of Dimensionality[END_REF]. Depending on the modeling framework, policies can be deterministic or stochastic with variable degrees of complexity. For many applications, however, a simplified framework can be formulated using "if-then" conditionals that instruct when to trigger an adaptation process [START_REF] Zhao | Decision rule-based method for flexible multi-facility capacity expansion problem[END_REF]. Under this formulation, the conditional threshold and the size of interventions become the main parameters to be defined. These elements can be used to construct more complex policies. For instance, a more complex policy can be constructed from the if-then conditional by defining the magnitude of the adaptation as the optimal change that minimizes the expected present value of the associated costs and maximizes the expected utility, while constraining the desired performance to remain at certain levels.

Generally, flexibility management requires not one but a complete sequence of decisions. The sequence starts with an initial decision to define the design characteristics of the system, i.e., the initial dimensions and the flexibility (maximum adaptation range and adaptation costs). Once the system is fully commissioned and starts operation, a monitoring process is required to verify the adaptation conditions defined by the policy. The successful implementation of any policy is conditioned on an inspection and monitoring program. When the conditions are met, a decision is made concerning the magnitude of the adaptation. The complexity of this problem lies in that every decision made will affect the future state of the system and, in consequence, the input for the policy in the future. Furthermore, the future state will also depend on uncertain external conditions that may not be stationary (changes in traffic, demand, climate change, etc.) [START_REF] Chester | Toward adaptive infrastructure: flexibility and agility in a non-stationarity age[END_REF]. Thus, the problem of managing flexibility is a classic example of sequential decision problems.

In summary, designing and managing flexible systems requires at least the following tasks (for one flexible component):

1. determine the optimal initial design x 0 , 2. determine the optimal initial flexibility level f 0 , 3. define the policy π to represent the management decisions, 4. following π, determine the optimal timing τ 1 , τ 2 , . . ., and size y τ1 , y τ2 , . . . of future adaptations.

The model presented in this paper addresses the first two items. Figure 2 presents a diagram that synthesizes this design and management process.

Sequential decision-making

Typical decision-making models require an agent and an environment. At time t, the agent observes the state of the environment s t and selects an action a t according to the policy π. Solving a sequential decision problem is not a unique problem; each model may have many different elements that need to be determined, which requires different methods. For instance, a key problem is finding the optimal policy, i.e., the optimal function that selects an action a t for each pair of action and state history a 1:t-1 , s 1:t-1 . Problems of optimal policy are usually solved using methods such as dynamic programming (DP), approximate dynamic programming (APD), reinforcement learning, genetic algorithms, and Monte Carlo tree searches. A different problem is finding the optimal initial decisions that optimize a certain performance measurement over the system lifetime, i.e., the optimal design that minimizes discounted adaptation costs. In this case, the optimal initial decisions are determined considering future sequences of decisions that may result from a particular decision rule and a probability model, even if the optimal sequence of decisions is not determined. Decision problems focused on finding a robust, optimal initial decision can be solved using MS programming [START_REF] Brandimarte | Numerical Methods in Finance in Economics[END_REF][START_REF] Shapiro | Lectures on Stochastic Programming[END_REF][START_REF] Birge | Introduction to Stochastic Programming[END_REF].

The problem of flexibility and the value of flexibility as presented before is better suited to be solved using MS programming. The next section presents an overview of the theory behind MSPs and 4 presents the proposed MS model for the problem of flexibility.

Stochastic Programming

Two-stage program with fixed recourse

The simplest stochastic programming model is the linear two-stages model with recourse shown in Equation 2 [START_REF] Birge | Introduction to Stochastic Programming[END_REF], where the objective is to minimize the costs at the initial stage t = 0 and at the recourse stage t = 1:

min x0 f = c T x 0 Initial stage costs + E ξ [min q(ω) T x 1 (ω)] Recourse costs s.t. Ax 0 = b H(ω)x 0 + Jx 1 (ω) = h(ω) x 0 ≥ 0, x 1 ≥ 0 (2)
At the initial stage, the decision x 0 ∈ R n1 is made considering that the first stage data is known: matrix A ∈ R m1×n1 and vector b ∈ R m1 . At the second stage, a random event ω ∈ Ω is realized (Ω represents the set of all random outcomes) and the second stage information is revealed:

technology matrix H(ω) ∈ R m2×n1 , recourse cost vector q(ω) ∈ R n2 , and right-hand side vector h(ω) ∈ R m2 . Then, the second stage decision x 1 (ω) ∈ R n2 , the recourse decision, is made. This model is classified as with fixed recourse because the recourse matrix J ∈ R m2×n2 is assumed deterministic. This assumption implies that the effect of the recourse decision x 1 in the constraints is known, which greatly simplifies the solution of the problem and still can be used to model a wide range of real life situations.

Each element in the random matrix and vectors H(ω), q(ω), and h(ω) is a random variable.

Therefore, the random vector ξ = (q(ω), h(ω), vec(H(ω))) with support Ξ ⊂ R d can be constructed to represent all the data that depends on the random event ω. By taking the expectation with respect to ξ, the estimation of the second stage costs is considering all possible realizations of ω (at least those with non-zero probability).

A key characteristic of stochastic programs (two-stages or multistage) is that the recourse decisions x 1 will be different for each realization of the random event ω. Therefore, if the random parameter distribution is discretized into r intervals, there will be r different x 1 solutions, which is not very useful from a planning perspective. The true value of stochastic programs lies in the initial stage solutions x 0 because they are unique for all trajectories of the random parameter. These decisions are being selected to guarantee that the expected value of future costs (which depend on ξ and x 1 ) will be minimal. The next section presents how this model can be extended to the multistage case.

Multistage stochastic programs

Multistage stochastic programming is a modeling framework that generalizes the two-stage problem to a sequence of recourse decisions. In the two-stage problem, the realization of the uncertain parameter becomes known at the second stage and the decision-maker has the option to make a second decision x 1 -the recourse-to adjust the initial decision x 0 . When the approach is extended to the multistage case, every recourse decision x t will be selected based on the expected value of this decision for stage t + 1 assuming that every future decision will be optimal [START_REF] Birge | Introduction to Stochastic Programming[END_REF]. This results in a particular sequence of recourse decisions for every sequence of uncertain parameters.

The general MS programming model can be formulated as follows [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF]:

min x1∈X1 f 1 (x 1 ) + E[ inf x2∈X2(x1,ξ 2 ) f 2 (x 2 , ξ 2 ) + E[. . . + E[ inf x T ∈X T (x T -1 ,ξ T ) f T (x T , ξ T )] . . .] (3)
where x t ∈ R nt is the vector of decisions variables at time t; ξ t ∈ R Mt is the vector of uncertain parameters at time t (with M t as large as n t + m t + n t × m t ); functions f t (x t , ξ t ) define the cost of making the decision x t ; and X t represent the set of constraints that define the feasibility regions.

By using an expectation functional in the recourse functions, it is being assumed that the decision-maker is risk-neutral. Risk-averse attitudes can be modeled by replacing the expectation functional with risk functionals that consider both the mean and the variance of the random outcome in the optimization process. Functionals such as mean-variance, semi-deviations, weighted mean deviations from quantiles, and average value at risk can be used to generate solutions with limited risk [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF]. Nonetheless, for the case of infrastructure networks where the performance is considered long-term and the structural integrity is not being subject to optimization, the expectation functional can be used as long as the number of scenarios considered is large enough for the Law of Large Numbers to apply. Other requirements defined by the decision-maker, e.g. performance levels, can be specified either as hard constraints in the feasibility set X t or the deviations can be penalized in the cost function.

In an MSP there may be four uncertain elements: future costs, right-hand side vector, technology matrix, and recourse matrix. The right-hand side vector is usually associated with demand, the technology matrix with the initial decision response, and the recourse matrix with the recourse decisions response. For instance, the technology matrix could represent the productivity of the initial conditions while the recourse matrix could represent the productivity of the recourse decisions, and both affect the system's ability to serve the random demand given by the right-hand side vector.

Clearly, the decision of which elements are considered random in the model depends on the type of system and the modeler's assumptions.

The decision on how to model the random parameters depends entirely on the type of problem.

Historical data, econometric projections, experts' opinions, all can be used to define the probability distributions of the elements selected to be random. The scenario generation techniques discussed in the following section are used to transform these probability distributions and the filtration structure of the information into discrete approximations in the form of scenario trees that can be used to solve numerically the MSP.

An important property of MSPs is that the information available is represented as σ-algebras

F t .
As the stages move forward, the information available increases, which is modeled as a sequence of increasing σ-algebras: F t ⊆ F t+1 . This growing sequence of σ-algebras is known as a filtration.

This is an important property of the problem because it implies that decisions at each stage x t can only consider the information available up to that stage [START_REF] Gröwe-Kuska | IEEE Bologna PowerTech -Conference Proceedings[END_REF]. This requirement can be explicitly or implicitly modeled in the constraint set and it is known as non-anticipativity constraints.

Challenges in stochastic programming

The MS programming approach, while useful to capture the presence of uncertainty in planning problems, faces limitations due to its computational complexity. The first limitation is defined by the representation of the uncertain parameters in the model. Regardless of whether the distribution of the random processes is known or if only sample paths are available, this information has to be transformed somehow to be entered as an input to the model. Scenario trees are the usual approach to represent the uncertain elements in MSP and constructing them implies the finite discretization of the original distributions (if they are continuous in the first place).

To properly represent the original distribution a large number of points are required. Furthermore, the nested structure of scenario trees results in exponential growth of the number of scenarios with the number of stages. These two issues (many points and exponential growth) render the problem intractable if a detailed representation of the underlying distribution is desired in a model with many stages. For this reason, methods for the generation and reduction of scenario trees constitute a very active area of research. Høyland and Wallace [START_REF] Høyland | Generating Scenario Trees for Multistage Decision Problems[END_REF] developed an approach based on minimizing a distance measure (e.g. squared norm) between a predefined statistical property and the statistical property from the approximation. This approach is known as moment matching and has been further expanded in various works (e.g. [START_REF] Xu | Scenario tree generation approaches using K-means and LP moment matching methods[END_REF]). However, as shown in [START_REF] Høyland | Generating Scenario Trees for Multistage Decision Problems[END_REF], matching some moments of the distributions is not enough to guarantee that the solutions of the stochastic program will be similar. For this reason, Heitsch and Römisch [START_REF] Heitsch | Scenario tree modeling for multistage stochastic programs[END_REF] developed scenario generation heuristics based on the concept of stability. Stability is a property in stochastic programs that guarantees that the optimal value does not change excessively by changing the scenario tree if both come from the same distribution. The approaches developed by Heitsch and Römisch [START_REF] Heitsch | Scenario tree modeling for multistage stochastic programs[END_REF] recursively reduce and bundle scenarios using the concept of filtration distance. The more recent approach developed by Pflug and Pichler [START_REF] Pflug | Dynamic generation of scenario trees[END_REF], based on the concept of nested distance, allows generating the structure of the tree dynamically to meet a prescribed precision.

The second limitation is associated with the solution of the optimization problem. The nested structure found in MSPs severely complicates the problem and special algorithms are sometimes required to provide a solution. Two general approaches are usually recognized in the literature [START_REF] Rosa | On augmented Lagrangian decomposition methods for multistage stochastic programs[END_REF][START_REF] Ruszczyński | Decomposition Methods[END_REF][START_REF] Bakker | A structuring review on multi-stage optimization under uncertainty: Aligning concepts from theory and practice[END_REF]: Primal and Dual decomposition. In the primal decomposition approach, the problem is divided into a master problem for the current stage and sub-problems for the previous stages.

During each iteration, cuts are generated to linearly approximate the recourse function and generate a candidate solution x * [START_REF] Birge | Introduction to Stochastic Programming[END_REF]. Examples of primal decomposition algorithms can be found in the works of Birge [START_REF] Birge | Decomposition and Partitioning Methods for Multistage Stochastic Linear Programs[END_REF], Ruszczyński [START_REF] Ruszczyński | A regularized decomposition method for minimizing a sum of polyhedral functions[END_REF], and Ruszczyński [START_REF] Ruszczyński | Regularized Decomposition of Stochastic Programs: Algorithmic Techniques and Numerical Results[END_REF]. In the dual decomposition approach, non-anticipativity constraints are relaxed and moved to the Lagrangian. The division between a master problem and many sub-problems also exists in this approach, but here the sub-problems are connected with scenarios instead of stages. As in the case of primal decomposition, most of the advances came during the late 80s and early 90s with the Progressive Hedging Algorithm [START_REF] Rockafellar | Scenarios and Policy Aggregation in Optimization Under Uncertainty[END_REF] and the Diagonal Quadratic Approximation [START_REF] Mulvey | A diagonal quadratic approximation method for linear multistage stochastic programming problems[END_REF]. Most recent efforts have focused on extending these approaches to mixed-integer problems, as in the Branch and Price approach from Lulli and Sen [START_REF] Lulli | A Branch-and-Price Algorithm for Multistage Stochastic Integer Programming with Application to Stochastic Batch-Sizing Problems[END_REF].

Nonetheless, not all the MSP required special algorithms to be solved. Alternative approaches exist to exploit particular structures of the problem. For instance, if Block-Separability Recourse (BSR) exists, the MSP can be analyzed as a two-stage problem with aggregate-level decisions being made at the first stage and planning-level decisions happening during the second stage. Commercial solvers are also capable of handling MSPs with linear constraints and non-linear objectives [START_REF] Dupačová | Scenarios for Multistage Stochastic Programs[END_REF] while state-of-the-art solvers can handle problems in the deterministic equivalent formulation [START_REF] Bakker | A structuring review on multi-stage optimization under uncertainty: Aligning concepts from theory and practice[END_REF].

Multistage Stochastic Programming Model for Flexibility

Problem description

An important challenge faced by infrastructure network systems consists of managing changes in demand such as traffic, consumption, or population growth efficiently. This can be accomplished with sporadic adaptations in size, capacity, or any other performance parameter to maintain minimum safety and operational standards. Decisions on the timing and sizes of these interventions should be carefully made to maintain an acceptable level of performance while keeping the costs minimal. In this section, an MSP has been adapted to model this problem.

Consider a network represented as a directed graph G(E, V ), where E is the set of edges and V a set of nodes, such that every edge e uv ∈ E connects nodes u and v with u, v ∈ V . Each edge is associated with a flow w t (u, v) occurring at time t. This flow can represent traffic, electricity, water, information, etc., moving through typical infrastructure network systems such as highways, power distribution networks, water supply and sewage distribution networks, and telecommunication networks. Simplifying the notation, each edge e can be characterized by a vector of design and operational variables x t,e ∈ R ne that may change over time (e.g., number of lanes in a highway or capacity in a power transmission network), where n e represents the number of variables considered

for edge e. These potential changes may occur as a result of some external random phenomena (e.g., change in traffic demand or flow), described by event ω ∈ Ω, and are represented by the variables y t,e ∈ R ne , in the case of flexible adaptations, and z t,e ∈ R ne for unplanned adaptations.

An additional element is the flexibility built into each design or operational variable, which is represented in the model as f 0,e ∈ R ne (a decision variable), the maximum range of change allowed by design for x t,e . This is not the flexibility vector fv presented in Equation 1, but the difference between x max and x 0 . For simplicity, f 0,e will be regarded as the flexibility from now on.

The flexibility decisions are made at the initial stage, together with the decision concerning initial design and operation properties x 0,e , and their purpose is to restrict the flexible adaptations y t,e that can occur in the future. These modeling decisions allow introducing the problem of balancing three approaches: investing at the initial stage to reduce future adaptation costs, investing at the initial stage to have robust/unchangeable systems, and not investing and risking larger future adaptation costs. This is, deciding the optimal values of x 0,e , f 0,e , y t,e , and z t,e . It is expected that most solutions will involve a combination of the three strategies. units within the flexible range f 0,e (i.e., the initial design value of variable i plus the sum of all the flexible adaptations over the lifetime must be less or equal than the initial value plus the flexibility x 0,e (i) + T t=1 y t,e (i) ≤ x 0,e (i) + f 0,e (i)), and d e (z t,e , ω) is the cost function of performing an adaptation z t,e outside the flexible range. This option is included in the model to represent the fact that any infrastructure network can be modified if enough resources are invested; however, this does not mean that the alternative will be efficient. Additional costs related with operation are the flow cost q(w e , ω), the operation and maintenance (O&M) cost g e (x t,e (i), ω), and the function of revenue received from exploiting the network h e (ω).

Definition of costs

Each decision variable

The parameters of these cost functions depend on, first, the type of function used for the representation and, second, the procedure performed to obtain the data. For many applications, a linear or convex function is enough to represent the required behavior. However, to represent specific conditions, e.g. economies of scale [START_REF] Cardin | Improving the Lifecycle Performance of Engineering Projects with Flexible Strategies: Example of On-Shore LNG Production Design[END_REF], concave functions may be needed. The values of these parameters can be obtained from autoregressive models, reports in the literature, experts' opinions, industry surveys, and sensitivity analyses. If the model considers random costs, then the cost function parameters become random variables and additional parameters for their distributions must be defined (using historical data, available literature, surveys, etc.).

In the model proposed, all the cost functions can be represented using convex functions. If a typical linear cost representation is used, then the only cost parameter is the unitary cost coefficient.

In this case, the only requirement is that the unitary costs of functions d e (unplanned adaptations) must be larger than the unitary costs of functions c e (planned adaptations); otherwise, the model will never assign value to the flexible option. If a different convex function is used, then it must be guaranteed that d e (x) > c e (x) ∀ x > 0. This does not mean that flexibility will always be preferred due to the additional cost be of introducing flexibility at the initial stage. However, if the condition is not met, then flexibility will never be a viable option from a cost standpoint.

MSP problem formulation

The elements described in the previous sections are combined in the following objective function b e (f 0,e (i))

Flexibility introd. cost

+ E ξ 1       min γ 1       e∈E ne i=1 c e (y
+ e∈E ne i=1 g e (x T,e (i), ω) ± h(ω) . . . ( 4
)
where i indexes the design and operational parameters being considered for each edge e, n e is the number of design and operational parameters, and γ t is the discount factor. To represent the limitations and special conditions associated with the management of flexibility and the optimal flow distribution inside the network, the following constraints are formulated:

First, the constraints for flow conservation between converging and diverging links at a node k are:

v w (k,v) - u w (u,k) = φ k (ω) ∀u, v, k ∈ V, (5)
where φ k = 0 if k is a transshipment node and φ k = 0 if k is either a supply or a demand node.

Second, the constraints for the minimum and maximum flow that can move through each link e are: where the lower and upper limits lower , upper are part of the design and operational variables vector x t,e and may have a flexibility value associated (can be adapted).

Third, the constraints that control the trigger of an adaptation process:

V (x t,e , ξ t ) ≥ π t,e t = 1 . . . T, ∀e ∈ E (7)
where the performance function V provides a quality measure of the design and operation variables x t,e in comparison with the pressure caused by the random parameters ξ t . This set of constraints establishes that the performance measure must be larger than a certain value π t,e given by the management policy. In this way, the constraints defines the criteria to make a change. For example, it can define the maximum ratio flow/capacity in an edge before it has to be expanded.

Fourth, the constraints for the maximum planned adaptation in a design or operation variable in an edge are:

t j=1 y j,e (i) ≤ f 0,e (i) t = 1, . . . , T, i = 1, . . . , n e , ∀e ∈ E (8)
These constraints establish that the sum of the total planned changes for a design/operation variable i in an edge e must be, at most, the flexibility built for that variable i and edge e.

Fifth, the constraints to update the state variables of an edge in accordance with the available actions are:

x t,e (i) = x t-1,e (i) + y t-1,e (i) + z t-1,e (i) t = 1, . . . , T, i = 1, . . . , n e , ∀e ∈ E (9)
Other constraints can be added depending on the nature of the network, or to properly characterize the adaptations (e.g. non-negativity constraints restrict adaptations to expansions; integrality changes the nature of the adaptations from modular to continuous). There are non-anticipativity constraints implicitly formulated in the other constraints to restrict the knowledge available to make a decision at each stage; this is, at each time t, the decisions can only depend on the information available up to t.

The decision stages t are selected considering the frequency of the decision process. This frequency may not coincide with the frequency of the random processes, in which case one decision stage may encompass many random process periods. For instance, it is not expected for an infrastructure network to be expanded every year, even if the demand is reported every month. In such a case, the decision stage t can consider the average, the maximum, or a quantile of the demand periods clustered in t. It is also not required for the decision stages t 1 , t 2 , . . . to be equal in size.

Considerations regarding the life-cycle of the network and the reliability of the uncertain data may result in a higher frequency of decisions during the first years of the network. It must also be considered that increasing the frequency of the decision stages also increases the size of the problem.

In brief, it is a modeling decision that depends on the purpose of the model and the experience of the modeler.

Approximate Solution of MSPs

Criteria for approximating the solution

Solving a multistage stochastic program using the original distribution of the uncertain parameters ξ t may not be feasible in many real applications because the set of possible paths (realizations of the demand) is infinite. For this reason, several methods have been developed to approximate the distribution of the random phenomena [START_REF] Heitsch | Scenario tree modeling for multistage stochastic programs[END_REF][START_REF] Gröwe-Kuska | IEEE Bologna PowerTech -Conference Proceedings[END_REF][START_REF] Löhndorf | An empirical analysis of scenario generation methods for stochastic optimization[END_REF]. Overall, the available methods focus on:

minimizing some measure of the distance between the original and the approximated distribution (e.g. Wasserstein distance) [START_REF] Heitsch | Scenario tree modeling for multistage stochastic programs[END_REF][START_REF] Pflug | Multistage Stochastic Optimization[END_REF], matching moments between the distributions [START_REF] Høyland | Generating Scenario Trees for Multistage Decision Problems[END_REF], and generating samples (e.g. Monte-Carlo sampling [START_REF] Shapiro | Monte Carlo Sampling Methods[END_REF]. As a result, the original distribution becomes a set of scenarios (each with an associated probability), usually organized in a tree structure. This scenario tree represents the increasing finite filtrations of available information, as shown in Figure 3.

Fig. 3. Scenario tree representation of increasing finite filtrations

The accuracy of the stochastic program solution highly depends on the quality of the resulting discretization [START_REF] Guastaroba | On the effectiveness of scenario generation techniques in single-period portfolio optimization[END_REF]. However, it is important to keep in mind that the main objective is to optimize the quality of the solution obtained by the stochastic program and not necessarily to produce an optimal discretization of the original distribution. This can be achieved by having stability and low bias in the solution [START_REF] Kaut | Evaluation of scenario-generation methods for stochastic programming[END_REF]. Stability means that the variation in the optimal value is minimal among different scenario trees (obtained from the same distribution). Low bias refers to a small gap between the real (unobtainable) and the approximated solution. The most direct approach to achieve both stability and low bias is to simulate a large number of scenarios [START_REF] Beraldi | Generating scenario trees: A parallel integrated simulationoptimization approach[END_REF][START_REF] Kaut | Evaluation of scenario-generation methods for stochastic programming[END_REF]. A trial and error procedure is usually necessary to achieve these objectives with the smaller scenario tree as possible.

Discretization method by clustering

As explained before, it is generally not possible to solve an MSP using the real distribution of ξ t ; instead, an approximation is required. Typical approximate distributions are represented in the form of scenario trees. In this paper, the discretization process is performed by combining Monte-Carlo sampling and a k-medoids clustering algorithm (see Figure 4).

The methodology requires listing the properties of the random process to be simulated, which can include values such as the mean and the standard deviation, or additional shape parameters.

Then, the time mission is divided in a finite set of time steps (stages) t = 1, . . . , T , and a large number N of random trajectories are generated (see Figure 4(a)). At the first time stage t = 1, the coefficient k 1 is defined to represent the number of clusters in which the trajectories will be aggregated; this coefficient is known as the branching factor. Then, the trajectories are combined into k 1 clusters using a k-medoids clustering algorithm. At the next time stage, the trajectories grouped inside each cluster are further divided into k 2 clusters using the same algorithm. This process is repeated until the end of the time mission when a total of K = T t=1 k t clusters are generated. This process results in a symmetric scenario tree.

This method requires a large number of trajectories per cluster to avoid clustering single trajectories. Considering that the number of scenarios grows exponentially with the number of stages, the number of required trajectories may be difficult to generate for large scale applications. While it works well for problems with a number of scenario paths of order 10 5 to be clustered in trees with a number of scenarios of order 10 3 , problems of larger scale may require more complex approaches. For instance, a larger number of clusters can be used at the initial stages to obtain a finer discretization due to the importance of the prediction at these early stages (Figure 4 number of variables and stages due to the limited number of characteristics that can be adapted and the low frequency of changes.

Approximated flexibility MSP model

The approximation of the distribution of ξ t using a scenario tree can be combined with the original flexibility MS approach ((4)) to produce an approximated MSP that can be solved with numerical solvers. By replacing the probability model with an approximate, simpler, model with finite support [START_REF] Pflug | Multistage Stochastic Optimization[END_REF], the original filtration becomes a finite filtration with a tree structure (scenario tree). This approximation leads to the formulation of an equivalent deterministic program, as shown in Equation 10:

+ e∈E q k (w T,k ) + e∈E ne i=1 g e,k (x T,e,k (i)) ± h T,k (10)
where s indicates the number of scenarios and p k is the individual probability of each scenario.

Scenarios are understood as a complete path from the root node to a leaf node in the scenario tree.

This formulation of the objective function clearly shows that a particular solution is obtained for each scenario k except for the initial stage variables (x 0 and f 0 ) which are the same for all scenarios.

Formulating the MSP in this form requires to explicitly develop the non-anticipativity constraints. The objective of these constraints is to ensure that scenarios with the same history up to stage t are indistinguishable [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF].

x t,e,j (i) = x t,e,k (i) y t,e,j (i) = y t,e,k (i)

z t,e,j (i) = z t,e,k (i) ∀ j, k for which ξ t,j = ξ t,k , t = 1, . . . , T, i = 1, . . . , n e , ∀e ∈ E (11)
To illustrate the non-anticipativity conditions, Figure 5 shows a fragment of a scenario tree. The figure shows that Scenario s 1,tr and Scenario s 2,tr are identical up to stage t q , same as Scenario s 3,tr and Scenario s 4,tr . Furthermore, the four scenarios are indistinguishable up to stage t p . The restrictions in Equation 11ensure that the decision variables comply with the non-anticipativity conditions given by the scenario tree of the random process.

The procedures described in Sections 5.1-5.3 to formulate and solve the proposed MSP are summarized in the chart presented in Figure 6.

Model scalability

The scalability of the proposed model can be analyzed for two elements: the scenario tree generation algorithm and the solution method. First, the scenario generation algorithm is based on the general approach described by Dupačová et al. [START_REF] Dupačová | Scenarios for Multistage Stochastic Programs[END_REF] where a set of scenario paths are generated according to a probabilistic model (based on, for instance, historical data series), a scenario tree structure is predefined according to some heuristic (e.g. detailed branching for early stages and coarse branching for later stages), and a clustering algorithm is applied to the scenario paths based on some dissimilarity measure. The scenario path generation process, even for the multivariate case, scales linearly in the worst case, and it is never regarded in the literature as a worrisome source of computational complexity. The clustering algorithms, however, do not scale as well due to the required comparisons between all the data points (distance or dissimilarity measures).

The k-medoid algorithms used in this paper are the Partitioning Around Medoids (PAM) and the Clustering LARge Applications (CLARA) implementations in Matlab ® . The CLARA algorithm is used to cluster the first 2-3 stages of the scenario tree where the number of data points is large (> 5000). This algorithm is of order O(k 3 ) [START_REF] Schubert | Faster k-Medoids Clustering: Improving the PAM, CLARA, and CLARANS Algorithms[END_REF] and can be used for large data sets to be bundled in Formulate the constraints that represent the adaptation policy (Eq. 7)

Implement the deterministic equivalent problem (Eqs. 5-11) using the parameters, cost functions, and scenario tree Solve the optimization problem and assess the quality of the solution using duality gap,VSS, etc. of points in an also small number of clusters, which can be done using the PAM algorithm of order

O(k(n -k) 2
). While these clustering algorithms work well for a wide range of data set sizes and a limited number of branching factors, the structure of the scenario tree requires that the algorithms have to be executed an increasing number of times at each stage. This is, the number of calls for the clustering algorithms grows exponentially with the number of stages. Therefore, the scalability of the scenario generation procedure depends on both the scalability of the clustering algorithms and the size of the scenario tree.

The second element that affects the scalability of the model is the solution method. As explained in Section 3.3, solving an MSP is challenging due to the considerable number of variables, the interdependence between them, and possible non-linear constraints. While special methods have been developed to solve MSPs, in some cases commercial solvers are enough to obtain a solution in reasonable times for problems with hundreds of thousands of variables. One of these cases is when the problem has a linear (or non-linear) objective and linear constraints [START_REF] Dupačová | Scenarios for Multistage Stochastic Programs[END_REF]. By introducing the scenario tree representation to formulate the deterministic equivalent problem, the result is a large scale linear (or non-linear) optimization problem that can be solved efficiently by commercial solvers. The model presented in Section 5.3 has linear constraints and the cost functions can be properly represented as linear or quadratic functions for many applications. Therefore, commercial solvers can be used to solve the proposed model, as is shown in the example in Section 6.

Nonetheless, the exponential growth in the number of variables as the number of stages increases poses a significant limitation to the frequency of the decision process. This limitation can be circumvented as the nature of infrastructure systems makes monthly or even annual adaptation plans unrealistic due to the time and resources required to complete an adaptation. The decision stages for infrastructure adaptation problems can encompass years and, in consequence, it is possible to model the complete or a considerable portion of the system's lifetime with a limited number of stages.

In addition to the linearity of the constraints and the control to the scenario tree explosion by limiting the number of stages, the proposed model has two additional properties that reduce its computational complexity. The first property is the diagonal structure in the constraint matrices.

Limiting the dependence of the decision variables to their immediate predecessor by adding state variables results in sparse matrices that are much faster to solve [START_REF] Brandimarte | Numerical Methods in Finance in Economics[END_REF]. The second property is block-separable recourse. This property allows transforming an MSP into a two-stage problem by separating the aggregate level decision variables from the planning level variables. The aggregate level decisions can represent capacity expansion decisions while the planning level decisions can refer to the use of this capacity [START_REF] Birge | Introduction to Stochastic Programming[END_REF]. If the aggregate level decisions do not depend on the planning level decisions, then the former can be moved to the first stage while the latter are placed on the second stage. Specifically, block-separable recourse exists if i) at stage t the cost functions can be written as:

f t (u t ) + g t (v t ) ( 12
)
where u t represents the aggregated level decisions, v t the planning level decisions, and f t and g t are the respective cost functions. And ii) if the constraints can be written as:

  T t 0 S t 0     u t-1 v t-1   +   W t 0 0 D t     u t v t   =   a t b t   (13)
where T t and S t are submatrices of the technology matrix, W t and D t are submatrices of the recourse matrix, and a t and b t are subvectors of the right-hand side vector. This equation clearly shows that both current aggregate and current planning level decisions depend only on previous aggregate level decisions.

The model proposed in Equations 4-11 is block-separable for the following reasons: i) the cost functions are clearly segregated between the aggregate level variables (x t , y t , z t , and f 0 ) and the planning variables w t ; ii) the constraints in Equations 5-9 and 11 are not enforcing any dependence between the set of aggregate decisions x t , y t , z t , and f 0 and the decisions w t-1 . Constraints 6 and 7 do generate a dependence between w t and x t-1 but this relationship is included in the definition of block-separability.

In summary, the proposed model has advantageous characteristics in its structure that can be exploited to vastly reduce its computational complexity and allow the use of commercial solvers to find a solution. The main bottleneck happens in the scenario tree generation procedure due to the combinatorial nature of the process. Despite these advantages, the number of variables in the problem still grow exponentially with the number of stages, and even if the MSP can be transformed into a large scale linear problem and a simple binary tree structure is used, having more than ∼ 22 stages implies hundreds of millions of variables.

Numerical Example

Case description

This section presents an example of modeling flexibility in network design and operation using the proposed MSP. The purpose of this example is threefold: i) to show how the model can be used to make decisions when managing flexible networks; ii) to study how the presence of flexibility affects the cost performance of the network; and iii) to identify the elements in the model that have the highest impact in the assessment of the value of flexibility.

The network used in the example is shown in Figure 7; it is a generic network that may describe telecommunication, energy, water distribution, or a transportation network. This example only considers one element in the random parameters vector ξ t = δ t to represent the demand that enters the network through node A and propagates throughout the network until it reaches node F , where it exits. This demand is modeled as a stochastic process defined by:

δ(t) = β 1 t + e β2t β 3 t sin(β 4 t) + B(t) (14)
where β 1 , β 2 , β 3 , β 4 are normally distributed random variables and B(t) is a Wiener process. Table 1 shows the values of the parameters used in the example. Figure 8 shows five realizations of this process.

β 4 N (0.9, 5.06 × 10 -2 ) A F D E C B x 0,(AB)
x 0,(AC)

x 0,(CE)

x 0,(BE)

x 0,(BD)

x 0,(EF)

x 0,(DF)

x 0,(DE) The design and operation variables vector x t,e of each link in the network only contains the link capacity x t,e . The initial capacity is notated as x 0,e , and together with the flexibility range f 0,e define the initial stage decisions. The maximum capacity that can be reached by flexibility for each link at any time is given by x max,e = x 0,e + f 0,e . However, the system can be modified beyond that value but at a larger cost. This may happen, for instance, due to larger times required to modify the system, larger material and labor costs that result from last-minute contracts, and the related disturbances and damages caused by the unexpected partial or total cease of operation. In the example, only capacity expansions are considered.

Input (demand)
The time mission of the network is T = 40 years starting from the date of commissioning. To reduce the size of the problem, it is assumed that adaptations can only occur at specific times: [START_REF] Roefs | Centralised, decentralised or hybrid sanitation systems? Economic evaluation under urban development uncertainty and phased expansion[END_REF][START_REF] Chester | Toward adaptive infrastructure: flexibility and agility in a non-stationarity age[END_REF][START_REF] Luss | Operations Research and Capacity Expansion Problems: A Survey[END_REF][START_REF] Singh | Dantzig-Wolfe Decomposition for Solving Multistage Stochastic Capacity-Planning Problems[END_REF][START_REF] Taghavi | A Multi-Stage Stochastic Programming Approach for Network Capacity Expansion with Multiple Sources of Capacity[END_REF][START_REF] Zhang | On the relation between flexibility analysis and robust optimization for linear systems[END_REF] (except for the sensitivity analysis in Section 7.1.2). This limitation allows downsizing the scenario tree, which highly impacts the computational performance. The reduction can be justified by the fact that infrastructure networks, and in general infrastructure systems, cannot be changed frequently in most of their dimensions without incurring high costs.

The model uses two families of cost functions. The initial design cost functions a e , flexible adaptation cost functions c e , and the unplanned adaptation cost functions d e (see Equation 10) are assumed linear of the form a e = α 0 x 0,e , c e = α f a y t,e , and d e = α ia z t,e , equal for all the edges. The linear cost assumption reduces the comparison to a straightforward relationship between the cost coefficients. Quadratic functions are used to represent the costs of introducing flexibility of the form b e = α f f 2 0,e . This non-linear increasing function is used to consider the effect of dis-economies of scale that may exist in systems that cannot be expanded indefinitely without restrictions, such as most of infrastructure networks.

Coefficient α 0 is equal to 5 cost units per capacity unit. Coefficients α f a and α ia vary for the analyses shown in Section 7.2. Operation and maintenance costs (O&M) g e is assumed equal to 2% of α 0 per unit of capacity installed, except for the last analysis. Flow costs q e are function of O&M costs g e as shown in Table 2.

Table 2. Flow costs summary

Link Flow Cost

A-B 0.05g e A-C 1.2q AB B-D 1.1q AB B-E q AB C-E 1.2q AB D-E 1.3q AB D-F 1.1q AB E-F 1.2q AB
The flow costs vary between links to represent an existing minimal cost route and increase the stability of the solution by avoiding a random allocation of initial capacity and flexibility.

Discretization of the demand space

The first task is the discretization of the demand space using the sampling and clustering algorithm presented in Section 5.2. This requires defining first the scenario tree branching factors which vary between [4 4 2 2 2 2] (256 scenarios) and [4 4 4 4 2 2] (1024 scenarios). Then, the algorithm generates between 15000 and 25000 demand trajectories according to Equation 14 depending on the total number of scenarios. The clustering algorithm processes these inputs to generate a scenario tree whose nodes contain the demand values at the decision stages [START_REF] Roefs | Centralised, decentralised or hybrid sanitation systems? Economic evaluation under urban development uncertainty and phased expansion[END_REF][START_REF] Chester | Toward adaptive infrastructure: flexibility and agility in a non-stationarity age[END_REF][START_REF] Luss | Operations Research and Capacity Expansion Problems: A Survey[END_REF][START_REF] Singh | Dantzig-Wolfe Decomposition for Solving Multistage Stochastic Capacity-Planning Problems[END_REF][START_REF] Taghavi | A Multi-Stage Stochastic Programming Approach for Network Capacity Expansion with Multiple Sources of Capacity[END_REF][START_REF] Zhang | On the relation between flexibility analysis and robust optimization for linear systems[END_REF] years, with an associated probability that depends on the number of trajectories clustered in each node.

MSP solution

To obtain the optimal network initial conditions, x * 0,e and f * 0,e for every edge e in the network, it is necessary to solve the multistage stochastic problem described by Equation 10, subject to the restrictions presented in Equations 5 to 11. The model assumes that every time a flexible adaptation y t,e occurs, the flexibility built into the system at the initial stage f 0,e is reduced until its depletion.

Once this point is reached only unplanned adaptations z t,(u,v) can be deployed. Finally, adaptations (i.e., increments incapacity) can only occur at t ∈ {0, 6, 12, 18, 24, 32, 40}.

The evolution of the network capacity is controlled by a policy π(ξ t , x t,e ) based on a safety criteria that establishes that the capacity of every edge at time t should be such that x t,e ≥ SF w t,e , with w t,e the flow moving through link e, and SF some safety factor. This formulation keeps the constraints in Equation 7 linear. Additional flow conservation constraints are formulated for each network node according to the topology shown in Figure 7.

The numerical results were obtained using the CPLEX 12.8 API for Matlab ® 2017b and the YALMIP toolbox as algebraic interpreter [START_REF] Löfberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF] with an Intel ® Core ™ i5-5200 2.20GHz processor. An average run for 1000 scenarios takes approximately 100 seconds of processing time, where 95% of the time is spent in the scenario tree generation algorithm and less than 5% in solving the optimization program.

The solver used (CPLEX) automatically selects the barrier optimizer due to the large size and sparsity of the problem. This algorithm stops iterating when the primal and dual solution are complementary, i.e., when a sum of products between the primal and dual solutions is smaller than a predefined tolerance. The tolerance of the complementary used in this numerical example is defined as 1 -8 .

Network Design for Flexibility: Results

Solving the MSP defined by Equations 5 to 11 a results in an optimal (minimal) value of the total present costs C * total , an optimal value of the initial capacity x * 0,e for each link e, an optimal value of the flexibility f * 0,e for each link e, and optimal values each variable y * t,e,k and z * t,e,k , that represent the evolution process of each network link, for every stage t and scenario k. Figure 9 shows three realizations of the demand process and the respective evolution of the whole network capacity.

Capacity/demand units Fig. 9. Network capacity evolution for three demand scenarios

This figure illustrates one of the main properties of MSPs: The optimal solution for the initial capacity, i.e. the initial stage decision, is the same for all scenarios while the recourse decisions that expand the system are particular for each scenario. By averaging the optimal costs of all the system's future responses, the model makes an initial decision that is optimal in average for all scenarios.

Considering that the solutions of the recourse decisions y * t,e,k and z * t,e,k are generated for each scenario, they are not suited for the comparisons intended in the forthcoming analyses. Only the initial solutions x * 0,e and f * 0,e are the same for all scenarios. For this reason, only these solutions, together with the optimal costs C * total , will be considered in the analysis of results.

Numerical stability of the model 7.1.1. Stability verification of the demand tree generation

Stability in the generation of the demand tree is essential for the accuracy of the model. If the program is not stable, the solution can vary considerably between scenario trees. To verify the stability of the solution, three scenario tree structures were generated using different branching factors. These factors determine how many nodes grow from each node at each time stage and, ultimately, the number of scenarios. Furthermore, the demand values in the nodes of the scenario tree vary due to the Monte-Carlo generation procedure. To address these sources of variability, the MSP was solved in batches of increasing size (from 2 solutions up to 40 solutions) and the resulting optimal costs were averaged for each batch. Figure 10 presents the coefficient of variation, CV , of the optimal costs for each batch size, and for three sets of branching factors (scenario tree structures). The scenario tree structures compared in Figure 10 were selected starting by the smallest possible symmetrical tree [2 2 2 2 2 2]. Additional nodes were increasingly added to the initial stages to improve the quality of the short-term predictions. The largest scenario tree considered, [4 4 4 4 2 2],

showed similar behavior with larger trees but with considerably lower computational times.

Figure 10 shows that there is variability due to the number of solutions being considered and due to the instability in the solution. As the size of the batch increases, the first source decreases and only instability remains. It is observed that the optimal costs are relatively stable with a CV = 2.5%. It is also observed that a larger number of scenarios reduces the number of solutions required to reach a stable CV (branching factor [4 4 4 4 2 2]). This result is expected as a larger number of scenarios represent better the underlying distribution, achieving higher stability as a result. Generally, only the optimal costs are compared in stability analyses and not the optimal solutions due to the typical flat objective functions found in MSPs [START_REF] Cario | Modeling and Generating Random Vectors with Arbitrary Marginal Distributions and Correlation Matrix[END_REF].

Planning horizon discretization into decision stages

The results from the previous section assumed that the planning horizon is divided into 6 decision stages. This section presents a sensitivity analysis to determine how this modeling decision affects the MSP results.

To this end, the planning horizon is divided into 4, 5, . . . , 10 decision stages, and the resulting MSP is solved. In the cases when the planning horizon is not a multiple of the number of stages, the initial stages encompass fewer years than later stages, but the purpose is to cover the entirety of the planning horizon as uniformly as possible. The scenario trees are designed to have the same number of scenarios (2 10 ) which requires using variable branching factors. Therefore, the discretization of the probability space for the intermediate stages will not be equal for all cases, with a difference of a factor of 2 at most in the number of scenario tree nodes for similar time instants. For instance, for 4 decision stages, there are 64 nodes representing year 20, while for 5 stages there are 32 nodes at year 16. Table 3 shows the optimal value, the optimal flexibility solution f * 0,AB for link AB, the optimal initial capacity solution x * 0,AB , and the processing time obtained for each number of decision stages.

The results in Table 3 show that when the number of stages that divide the planning horizon is 5 or larger, only the processing times are affected. For 4 decision stages, the model solutions seem to be significantly different, beyond what could be explained by the random variation present in the solutions for the remaining number of stages. This tendency was found in multiple executions. This difference could be explained by the relationship between the frequency of the required adaptations and the frequency of the decision stages. As shown in Figure 9, the time between adaptations can vary between 8 and 12 years on average, with longer times during the first half of the planning horizon. This means that during the system's lifetime only 3 or 4 adaptations are expected. Therefore, when the planning horizon has 4 decision stages or fewer, it can be expected that an adaptation will take place at every stage and that the magnitude of these adaptations (and the initial capacity) will be larger to accommodate the changes in the demand during such long time periods.

In summary, the discretization of the planning horizon into decision stages can affect the solution when there is insufficient opportunity to use the recourse to respond to the changes in the external conditions.

Out-of-sample behavior of the model

The verification of the out-of-sample behavior consists in solving the optimization program with a particular tree and using the optimal solution, i.e., x * 0,e and f * 0,e , as the initial solution for MSPs with different scenario trees of the same size and generated from the same stochastic process. Figure 11 compares the optimal costs obtained in the baseline program (horizontal line), with the optimal costs found in 10 out-of-sample programs. This figure shows that 70% of the out-of-sample costs have a difference of 1.2% or less with the base case, while the remaining 30% has, at most, a difference of 3.3%. The results from the stability analysis and the out-of-sample analysis suggest that even if there is some instability because the scenario tree is not infinite, it is under control and the program exhibits both in and out sample stability, which is a requirement for low bias and general stability.

Value of the stochastic solution

A second key element in the analysis of MSPs is the value of the stochastic solution (VSS).

This metric compares how much is gained by solving the stochastic program instead of solving a deterministic program where the scenario tree is replaced with the expected value of the random process. To calculate the VSS, first, the expected value problem is solved and the optimal solution is used as the initial solution for an MSP. The optimal costs obtained for this new program are known as the expected value solution (EV). The VSS is obtained as the difference between the original MSP optimal costs and the EV.

Using the scenario tree structure regarded as the most stable from the analysis in Section 7.1.1, [4 4 4 4 2 2], generating 50 scenario trees and averaging the VSS, it was obtained a value that is 16% of the average optimal costs. This result is not negligible and suggests that the VSS is large enough to justify solving a stochastic program instead of the much simpler expected value problem.

The value of flexibility

The third element of the analysis focuses on the effects of the cost functions on the value of flexibility. Different relationships between α 0 , α f a , and α ia are considered to determine the effect on: i) the optimal initial capacity x * 0,e , ii) the optimal flexibility f * 0,e , and iii) the optimal costs.

Figures 12 and13 show, respectively, the evolution of the optimal solution for the capacity x * 0,e and the flexibility f * 0,e in each link for different ratios α f a /α 0 . The critical route is highlighted with the dotted lines and light shadow. As expected, the larger the cost of future adaptations α f a compared with the cost of building the required capacity from the beginning α 0 , the less flexibility is introduced, and the system relies more on the initial capacity to serve the demand. This means that reducing α 0 or increasing α f a decreases the value associated with flexibility, and the system manager is willing to pay less to have it in the network.

It is interesting to observe in Figures 12 and13 that the optimal solution for the initial capacity

x * 0,e results in larger values for the links in the minimum cost route, while the optimal solution for flexibility f * 0,e is similar for all the links in the network, excluding DE. This suggests that the short and medium-term demand is mostly met by the minimum cost route until a point where it is less expensive to use the flexibility built in other links to expand them and reroute part of the demand.

Besides, the quadratic cost function used for the introduction of flexibility discourages the approach of adding too much flexibility into one link. is not as cost-effective as just relying on the initial capacity for the short-term and unplanned adaptations for the medium and long-term. This could explain the monotonic decreasing tendency for flexibility while the initial capacity behaves asymptotically. This means that reducing the cost of unplanned adaptations can only reduce the dependence on initial capacity up to the point where there is enough capacity to attend the short-term demand that could not be served in any other way. This point is dictated by the initial demand and its growth during the first decision period.

A final analysis was conducted to examine the effect of the unitary O&M costs on the optimal initial capacity and flexibility. By increasing the unitary O&M costs from 2% of α 0 to 10% of α 0 , the optimal costs increase approximately 20%. The optimal initial capacity in all links is reduced by approximately 12% which is compensated by a similar increment in the optimal installed flexibility.

If the O&M costs are further increased to 16% of α 0 , the optimal initial capacity decreases another 9% with a similar increment in the optimal flexibility. Increasing the O&M costs discourage large system configurations in the short and medium-term, which makes flexibility more valuable. By modeling the flexibility built at the initial stage as a decision variable that has an associate cost and that limits future planned adaptations, the program can determine the optimal configuration that minimizes costs. This formulation also provides an indirect measure of the value of flexibility by finding the optimal amount of flexibility that is cost-effective to introduce.

An approximate solution method was implemented using Monte-Carlo simulation and k-medoids clustering to generate the discrete representation of the stochastic process (scenario tree). The deterministic equivalent program was formulated and tested for a generic network with one input and one output node for a stochastic demand. The numerical solution showed stability for reasonably sized scenario trees and a not negligible VSS of around 16%. Furthermore, the results showed how the value of flexibility is affected by the discount rate and the ratios between all the costs involved, not just the cost of deploying planned adaptations. If the ratios α f a /α 0 , and α ia /α 0 are too large, the system will prefer to meet the short and medium-term demand with the initial capacity, lowering the value of flexibility. As these ratios increase, the adaptations are further pushed into the future, especially for low discount rates. In the case where the ratio α f a /α ia is large, unplanned adaptations become cost-effective, also decreasing the value of flexibility. In contrast, for a given α f a , as the ratio α ia /α 0 grows, flexibility becomes more desirable. These examples show how the value of flexibility constantly evolves due to complex decisions that happen at different time instants.

The proposed model and the procedure implemented in the numerical example assumed that all decision and state variables are continuous. A mixed-integer formulation may require specific algorithms (e.g. primal decomposition) to find a solution efficiently. Furthermore, all cost functions were assumed convex, which is a reasonable assumption for most applications, but if concave functions are required, the solution procedure may not efficient nor may guarantee to find a solution.

The scenario generation procedure can be extended to the multivariate case (for independent variables) but the combinatorial nature of the problem restricts the number of decision stages than can be modeled. The formulation and solution of the MSP is also restricted by the exponential scalability. The expectation functional used in the model (Equation 4) enforces a risk-neutral approach that may not be appropriate for some applications. Further works should consider a formulation of the model with risk-averse functionals.

Acknowledgements

The financial support from the Vice-rectory for Research at Universidad de los Andes and from

Fig. 1 .

 1 Fig. 1. Flexibility vector representation

  of the model has a cost function associated. These cost functions can represent costs related with activities of design, construction, adaptation, operation, and maintenance of the network edges. The model considers the initial stage costs (design and construction) as known, while the costs associated with future activities (adaptation, operation, and maintenance) can be affected by uncertainty. This distinction is represented in the following cost functions: for the initial stage, a e is the cost function of building edge e with the initial design parameters x 0,e , and b e is the cost function of adding features to network edge e that facilitate future changes, i.e., the cost of adding flexibility to the system. For example, b e can represent the additional construction costs for a larger foundation in a building to facilitate future expansions. For the recourse stages, c e (y t,e , ω) is the cost function of adapting the network edge parameters by y t,e

Fig. 4 .

 4 Fig. 4. (a) realization of system path trajectories; (b) clustering of sample paths at different stages.

  k (x 1,e,k (i)) ± h 1,k + . . . k (y T,e,k (i)) + e∈E ne i=1 d e,k (z T,e,k (i))

Fig. 5 .

 5 Fig. 5. Non-anticipativity conditions in scenario tree

  paths in the nodes of the scenario tree using the k-medoids or any other clustering algorithm Calculate the probabilities of each scenario in the tree from the number of paths clustered in each leaf node Formulate the constraints that depend on the network structure (Eq. 5 and 6)

Fig. 6 .

 6 Fig. 6. General procedure to formulate and solve proposed MSP model

Fig. 7 .

 7 Fig. 7. Network used in the numerical example.

Fig. 8 .

 8 Fig. 8. Realizations of the stochastic process

Fig. 10 .

 10 Fig. 10. Solution variability as function of number of scenarios and runs

Fig. 11 .

 11 Fig. 11. Out-of-sample solution comparison

α f a /α 0 = 2 Fig. 12 .Figure 14 2 Fig. 13 .Fig. 14 .Fig. 15 .

 212142131415 Fig. 12. Initial capacity evolution for different ratios α f a /α 0

Fig. 16 . 8 .

 168 Fig. 16. Installed capacity and flexibility in link AB as function of ratios α f a /α ia , α 0 /α ia and discount rate

the

  Regional Council of 'Pays de la Loire' within the framework of the BUENO 2018-2021 research program (Durable Concrete for Offshore Wind Turbines) is gratefully acknowledged. List of symbols and abbreviations a e Cost function of building edge e with intial parameters x 0,e α 0 Unitary cost of building initial design x 0,e α f a Unitary cost of performing flexible adaptation y t,e α ia Unitary cost of performing unplanned adaptation z t,e b e Cost function of adding f 0,e flexibility to edge e c e Cost function of changing edge e by y t,e CV Coefficient of variation d e Cost function of changing edge e by z t,e γ t Discount factor at time t G Directed graph dr Discount rate e Edge/link in directed graph G f 0,e Flexibility range of design/operation parameters for edge e w t,e Flow thought edge e at time t fv Flexibility vector g e Operation and maintenance cost function for edge e h(ξ t ) Revenue function MDP Markov decision process MS multistage stochastic MSP multistage stochastic program O&M Operation and Maintenance p k Proability of scenario k π Policy function q e Flow cost function for edge e ξ t Vector of random, external parameters at time t ROA real options analysis VSS Value of stochastic solution x t,e Vector of design/operation parameters at time t for edge e x 0,e Vector of initial design/operation parameters for edge e y t,e Vector of flexible adaptations at time t for edge e z t,e Vector of unplanned adaptations at time t for edge e 10. References

  1,e (i), ω)

								ne
						+		d e (z 1,e (i), ω)	+	q(w 1,e , ω)
							e∈E	i=1	e∈E
				Flexible adaptations cost	Unplanned adaptations cost	Flow cost
								
	+	e∈E	ne i=1	g e (x 1,e (i), ω)	±	Exploitation cost/revenue  h(ω)    	+ . . .
			O&M cost

+ E ξ T min γ T e∈E ne i=1 c e (y T,e (i), ω) + e∈E ne i=1 d e (z T,e (i), ω) + e∈E q(w T,e , ω)

Table 1 .

 1 Random demand process parameters

	Parameter	Distribution(µ, σ 2 )
	β 1	N (450, 1.26 × 10 4 )
	β 2	N (-0.02, 2.5 × 10 -5 )
	β 3	N (100, 625)

Table 3 .

 3 Sensitivity analysis results

	Decision stages Optimal value f * 0,AB	x * 0,AB	Processing time(s)
	4	361619	2546 14787	45
	5	340954	2657 13341	62
	6	334664	2614 12798	85
	7	339480	2661 12894	95
	8	329768	2793 11175	92
	9	341655	2664 12618	90
	10	341035	2693 12961	110

Effect of the discount rate

A similar analysis was conducted to examine the effect of the discount rate on the optimal values x * 0,e and f * 0,e . Figure 15 compares the capacity and flexibility of link AB (other links follow a similar tendency) for three discount rates: dr = {3%, 5%, 8%} as a function of α f a /α 0 . While the values of flexibility are mostly unaffected by the change in the discount rate dr, the initial capacity presents a clear decreasing tendency with an increment in the rate. In this case, larger discount rates make future adaptations (inside or beyond the flexible range) more attractive in the medium and long-term.

The decision of adding flexibility depends both on the cost of adding flexibility, which is incurred at the initial stage and it is not affected by the discount rate, and the costs of performing adaptations, which are discounted. Increasing the discount rate reduces the present value of all future cost making both planned and unplanned adaptations more attractive, or at least moves forward the time stage where they are the optimal option. As a result, the short and medium-term demand is served using planned adaptations (flexibility) and long-term demand is served with unplanned adaptations, in detriment of the initial capacity. Therefore, it is expected that for typical ranges of the discount rate (∼ 3 -10%) flexibility will be considered valuable. However, it could be argued that for unrealistically large discount rates the model would rely entirely on unplanned adaptations, except for the initial demand, abandoning completely the option of flexibility.

Costs of adaptation beyond the flexible range

In this analysis, the effect of the cost of unplanned adaptations α ia is determined by using different ratios α f a /α ia and α 0 /α ia to solve the MSP. Figure 16 shows a similar tendency as Figure 15, with the optimal initial capacity and flexibility decreasing as the cost of unplanned adaptations α ia decreases, and the initial capacity being affected by the increasing discount rate, which further bolsters the cost-effectiveness of unplanned adaptations. However, in this case, a reversing tendency is detected for the 8% discount rate. Upon closer inspection, it is revealed that there is a point where the short-term demand must be met before unplanned adaptation can take place. At this point it does not matter how much further α ia can be reduced; there is a minimum demand that must be met.

For lower discount rates it was cost-effective to use more flexibility to meet this short-term demand; however, for large discount rates, the initial cost of flexibility plus the adaptation cost