Title Hydrological Scenario Reduction for Stochastic Optimization in Hydrothermal Power Systems Authors Ignacio Aravena and Esteban Gil Universidad T´ecnica Federico Santa Mar´ıa, Valpara´ıso, Chile Journal: Applied Stochastic Models in Business and Industry, vol. 31, no. 2, pp. 231-240, March/April 2015. URL: http://onlinelibrary.wiley.com/doi/10.1002/asmb.2027/abstract DOI: 10.1002/asmb.2027. Abstract Most of the methods developed for hydrothermal power system planning are based on scenario-based stochastic programming and therefore represent the stochastic hydro variable (water inflows) as a finite set of hydrological scenarios. As the level of detail in the models grows and the associated optimization problems become more complex, the need to reduce the number of scenarios without distorting the nature of the stochastic variable is arising. In this paper we propose a scenario reduction method for discrete multivariate distributions based on transforming the moment-matching technique into a combinatorial optimization problem. The method is applied to hydro inflow data from the Chilean Central Interconnected System (SIC) and is benchmarked against results for the optimal operation of the SIC determined with the selected subsets and the complete set of historical hydrological scenarios. Simulation results show that the proposed scenario-reduction method could adequately approximate the probability distribution of the objective function of the operational planning problem. Keywords Data analysis, Scenario reduction, Stochastic modeling, Energy resources, Hydrothermal power generation, Mean square error methods, Moment methods (c) 2014 John Wiley & Sons, Ltd. Personal use of this material is permitted. Permission from John Wiley & Sons, Ltd. must be obtained for all other uses.

1

This is the accepted version of the following article: Aravena, I. and Gil, E. (2014), Hydrological scenario reduction for stochastic optimization in hydrothermal power systems. Appl. Stochastic Models Bus. Ind.. doi: 10.1002/asmb.2027, which has been published in final form at http://onlinelibrary.wiley.com/doi/10.1002/asmb.2027/abstract.

Hydrological Scenario Reduction for Stochastic Optimization in Hydrothermal Power Systems Ignacio Aravena Esteban Gil Universidad T´ecnica Federico Santa Mar´ıa, Valpara´ıso, Chile

Abstract Most of the methods developed for hydrothermal power system planning are based on scenario-based stochastic programming and therefore represent the stochastic hydro variable (water inflows) as a finite set of hydrological scenarios. As the level of detail in the models grows and the associated optimization problems become more complex, the need to reduce the number of scenarios without distorting the nature of the stochastic variable is arising. In this paper we propose a scenario reduction method for discrete multivariate distributions based on transforming the moment-matching technique into a combinatorial optimization problem. The method is applied to hydro inflow data from the Chilean Central Interconnected System (SIC) and is benchmarked against results for the optimal operation of the SIC determined with the selected subsets and the complete set of historical hydrological scenarios. Simulation results show that the proposed scenario-reduction method could adequately approximate the probability distribution of the objective function of the operational planning problem.

Keywords: Data analysis, Scenario reduction, Stochastic modeling, Energy resources, Hydrothermal power generation, Mean square error methods, Moment methods

1

Introduction

A major goal in stochastic optimization models of hydrothermal power systems is to address the water inflow variability and the risks associated to drought years. Therefore, optimization models need to be matched to an adequate representation of the hydrological variability so that the hydro uncertainty is properly represented while keeping a manageable problem size. DOI: 10.1002/asmb.2027

ASMB 2014

The operation of the Chilean Central Interconnected System (Sistema Interconectado Central, SIC) strongly depends on seasonal weather conditions, which causes that the hydro energy injected in the system varies from 40% to 70% of the total generation, between dry and wet years. There are 48 years of historic hydrological data [1] for the relevant inflows. However, using the entire dataset for some of the larger stochastic optimization problems (such as capacity expansion planning problems [2, 3]) is usually impractical due to the size and computational complexity of the mathematical problem [4, 5]. In this context, although there is a trade off between accuracy of the representation of the stochastic variable and the computational size of the problem, improvements on scenario reduction techniques can be very useful to consider the stochastic nature of the optimization problem while controling its size. Furthermore, as the level of detail in the models grows and the associated optimization problems become more complex, the need to reduce the number of scenarios arises. This reduction of scenarios should not distort the nature of the stochastic variable. In power systems, most applications of scenario reduction methods focus on unit commitment and short-term operation [6]. However, there have been some developments in using them for operational planning in hydrothermal power systems. For instance, in [7] the authors proposed a method for predicting and generating reduced scenarios trees for water inflows based on the Principal Component Decomposition (PCD) technique along with Periodic Autoregressive (PAR) models. Generally speaking, several different techniques to perform scenario reduction exist once the probability distribution of the stochastic variable is already known or approximated [8] (usually by a large number of scenarios). A commonly used approach is clustering the scenarios with the K-means algorithm, whereas other techniques specifically designed for stochastic programming use probabilistic metrics to ensure bounded differences in the optimal values of the decision variables between the original and the approximate probabilistic distribution [9, 10]. A less used approach consists of the moment-matching technique, aimed to approximate the probability distribution by matching its statistical properties [11]. Recently, [12] presented methods based on K-means and moment-matching to generate scenarios trees directly from the data, avoiding the creation of a large number of scenarios to later reduce them. In the Chilean SIC case, regulation requires the use of historical hydro data to represent hydroelectric generation variability. Thus, it is convenient that any simplification of the original hydrological dataset (scenario reduction) preserves at least a few of the original hydrological series for 2

comparison purposes. This goal can be accomplished using a Sub-Set of the Historical Scenarios (henceforth SSoHSc) to represent the Complete Set of Historical Scenarios (henceforth CSoHSc). However, this approach implies losing a degree of freedom in the scenario reduction, since the elements of the SSoHSc have to be selected from the elements of the CSoHSc instead of creating them anew to match the statistical properties of the CSoHSc. In this paper, we propose a methodology to select a SSoHSc (and their respective weights) for operational planning and capacity expansion planning of the SIC. Although in this study SIC data was used, the method can also be applied to other systems with similar characteristics. The main idea behind the proposed method is to turn the moment-matching technique into a combinatorial optimization problem, so the differences between the moments of the CSoHSc and the SSoHSc are minimized. The quality of the selection is evaluated a priori using the moment differences in the multivariate energy inflows distribution and a posteriori with one year of operation results obtained for the CSoHSc. The operational results for the CSoHSc are used to benchmark the SSoHSc in terms of moment differences for the generation cost. This paper is structured as follows. First, section 2 briefly describes an analysis of the CSoHSc that supports its representation using a SSoHSc. Section 3 gives a brief description of the moment-matching scenario reduction method and presents its adaptation for the selection of the SSoHSc. Section 4 presents the results of the selection for subsets of different size and shows how the moments of the SSoHSc compare to the CSoHSc. Section 5 applies the selected subsets of different size to an operational model of the SIC in order to assess how well the SSoHSc will perform in a concrete application, by evaluating how much the approximation of the original probability distribution made by the SSoHSc will distort results in the optimization. Finally, section 6 summarizes the paper’s main conclusions and discusses possible improvements and applications of the methodology.

2

Analysis of the historic hydrological dataset (CSoHSc)

The dataset [1] consists of series of monthly averages of the water inflows for each run-of-river (RoR) generator and weekly averages for each generator with a reservoir, from april-1960 to march-2009. This series can be converted (using the respective generators’ efficiencies) into total monthly/yearly energy inflows. As the energy inflows are variables directly related with system operation costs, they are more convenient to use for the present analysis (and

3

further subset selection) than the water inflows. Although the total monthly energy inflow series exhibits a yearly periodicity, it seems stationary in the long term. It presents strong autocorrelations with the two previous months and with the same month in the previous year. The yearly energy inflow series shows a large variance (standard deviation is around 19% of the average), a strong negative skewness due to extreme low values and a moderate kurtosis (mesokurtic shape). There are not any significant autocorrelations in the yearly energy, meaning that regardless whether a year is dry or wet it will have no effect on the following years’ hydrology. The analysis was conducted using R [13], and allowed us to conclude that: (1) The intra-annual (monthly) autocorrelations can be represented separately for each year and (2) since the yearly inflow energy appears to be stationary in the long-term and do not show any significant autocorrelations, there is no need of using a scenario tree or series of hydrological years with dependence. Therefore, the CSoHSc can be represented by a SSoHSc, as long as the average, variance, skewness and kurtosis of the full collection can be maintained (approximately) in a subset of years by selecting them properly, and by choosing a convenient weight for each year in the subset. It is important to consider at least up to the fourth moment to represent the negative skewness and the presence of extreme values in the full collection. It is also important to consider the autocorrelations and cross-correlations between the inflows. If the inflow variable is conceived as a vector X, in which each component is the total energy inflow in a month of the year, then the covariance matrix represents the autocorrelations of the process. Additionally, if the energy inflow is split into the respective inflows, the covariance matrix represents both the auto and cross-correlations between the inflows of the system. Therefore the covariance matrix should be considered for the selection of the years and their respective weights.

3 3.1

Adaptation of the moment-matching method Overview

The moment-matching method was first proposed by Miller and Rice in 1979 as a way to approximate univariate continuous distributions with discrete distributions. Except for the extension to multivariate distributions and the use of transformations in intermediate steps, the method remains the same. In the univariate case, the objective of this method is to build a discrete approximation (i.e. a discrete-valued variable associated with discrete 4

probability values or weights) that has exactly the same moments than the original distribution. This is accomplished by stating a system of equations with the moments of the original distribution on the right-hand-side and the expressions of the discrete distributions moments on the left-hand-side, in terms of the discrete values and their probabilities of occurrence. The number of possible values of the discrete variable determines the number of moments that can be met and the number of equations to state. The solution of the system of equations provides the discrete approximation of the original distribution. Before discussing how we adapted the moment-matching method to the problem, it is necessary to define some operations and variables. First let the component power operation of a vector (·)·i be defined according to the expression (1). Then the zero centered moments of order i of an equally likely multivariate discrete distribution with nJ = |J| possible values1 xj can be determined with the expression (2) and the moment of a subset Jr (Jr ⊂ J, Jr 6= J) associated with the probabilities ρj is given by (3).    i a1 a1  a2   ai     2 a =  .  → a·i =  .  (1)  ..   ..  aiN

aN µ(i) = m(i) =

1 X ·i xj nJ j∈J X ρj x·ij

(2) (3)

j∈Jr

The moments defined by (2) and (3) contain information about the probabilistic distribution of each component of X, but none about the correlation between pairs (or tuples) of components. In order to take into account the information of the cross-correlations we use the matrix of bivariate moments M (analogous to the elements in covariance matrix, but centered on zero instead of the respective averages), defined in (4). Its components are determined via the equations (5) and (6), where xst corresponds to the tth component of the xs vector. 1

The cardinality of the set A is denoted by |A|.

5



µ11 µ12  µ21 µ22   .. ..  . . M=  µk1 µk2   .. ..  . . µN 1 µ N 2 µkp = mkp =

··· ··· .. .

µ1p µ2p .. .

··· ···

···

µkp .. .

···

µN p

··· .. . ···

1 X xjk xjp nJ j∈J X ρj xjk xjp

 µ1N µ2N   ..  .   µkN   ..  .  µN N

(4)

(5) (6)

j∈Jr

Back to the problem formulation, if it is assumed that the possible values for the stochastic vector X are all linearly independent (which indicates that X is on Rk , with k ≥ nJ ; a common situation when working with real world data), then to match even the first moment is impossible because µ(1) is a linear combination of the nJ values of the original distribution, while the expression for m(1) contains only a subset of the nJ values and therefore µ(1) is out of the range of possible values for m(1) (µ(1) is out of the range space of {xj | j ∈ Jr }). A similar reasoning can be applied to higher-order moments. Since in the multivariate case it is impossible to perfectly match the moments of the original distribution, it is desirable that the moments of the subset be at least similar to the original ones. The problem of selection then becomes an optimization problem that can be expressed as (7), where the γ factors are used to rescale the difference of the moments into the same range and to properly weigh each difference (differences on the first moment are much more important than differences in the fourth moment) and K is the set of the indexes of pairs of components, i.e. K = {1, 2, . . . , N } × {1, 2, . . . , N }, N = dim(xj ).

6

min

Jr ,ρj

s.t.

4 X i=1

2

γi µ(i) − m(i) + γˆ2

X

2

X (k, p) ∈ K k>p

(µkp − mkp )2 (7)

ρj = 1

j∈Jr

ρmin ≤ ρj ≤ 1, ∀j ∈ Jr

Jr ⊂ J

The second term in the objective function corresponds to the difference in the bivariate moments, which are specified by pairs of components instead of a matrix since M is symmetric and, therefore, the elements in the triangular sub-diagonal (or super-diagonal) matrix of M are the only ones of interest for the objective function. The differences in the diagonal elements of M are already considered in the first term of the objective function. Assuming that the norms of the moment differences are related similarly to the norms of the moments, the weights of the objective function are calculated according to (8) and (9), where the factor σ was determined by evaluating the parts of the objective function for some test subsets giving adequate weights for values in the range 2.2 − 2.4. By adequate weights we mean that they cause the differences of the moments multiplied by the weights γi to be on the same order of magnitude. Thus, the participations in the objective function were around 60%, 20%, 10% and 10% for the first, second and bivariate, third, and fourth moments, respectively.

γi = γˆ2 =

(1) σ−1

µ

1

µ(i) σ 1 γ2 N −1

(8) (9)

The aim of the relation (9) is that the sum of the squared errors in the bivariate moments has a similar importance than that of the objective function (7) compared to the sum of the squared errors in the univariate second order moment. Finally, a set of constraints force the sum of the weights to be 1 and that each of the weights be at least ρmin in order to avoid selecting scenarios with very small probabilities of occurrence.

7

3.2

Implementation

The resulting problem (7) is a non-linear combinatorial optimization problem. There are several approaches for solving this kind of problems, most of them in the domain of heuristics and meta-heuristics. Depending on the problem structure there are a few exact solution algorithms. Algorithms designed specifically for minimizing the sum-of-squares in clustering problems [14] are of particular interest as they are structurally very similar to (7). Although that problem is simpler than the problem stated in (7), it can be used as a reference in terms of solution time and complexity, or it could even be modified to solve (7). In this study the problem (7) is solved by decomposing it into two nested problems; it can be noted that once the subset Jr is fixed, the non-linear combinatorial optimization problem (7) turns into the quadratic programming problem (10) that has to be solved for the weights ρj , ∀j ∈ Jr ; therefore the problem can be solved using the following scheme: X X X ajj · ρ2j + min ajl · ρj ρl + bj · ρj + c ρj

s.t.

j∈Jr

X

j, l ∈ Jr j
j∈Jr

(10)

ρj = 1

j∈Jr

ρmin ≤ ρj ≤ 1, ∀j ∈ Jr

Jr ⊂ J

• At the top level a meta-heuristic technique (Differential Evolution [15]) selects a test subset of years. • For each test subset the weights are calculated and the minimum value of the objective function is found using spectral methods [16]. The overall process is depicted in Fig. 1. As the evaluation of the test subsets can be very time consuming, before starting the meta-heuristic search the annual inflow information is partially aggregated, selecting the components to be aggregated by PAM2 [18]. The PAM algorithm is also used to pre-group the elements in the full scenario collection for both providing some smoothness to the subset Jr over J [8] and reducing the search space. 2

Ideally this aggregation should be carried out using Principal Component Decomposition (PCD) [17] but since in the case of study there are more components than elements, the PCD results indefinite.

8

The meta-heuristic search, specifically the evaluation of the subsets, is performed in parallel over several cores or processors, allowing this algorithm to take advantage of multi-core and multi-processor architectures. Using parallel computing allowed the reduction of simulation times by up to 85%. The exercise of selecting test subsets and calculating the minimum value of the objective function is repeated until a stopping criterion is met. For this case study, the stopping criterion was defined as 100 iterations of the metaheuristic search without improving the objective function (each iteration corresponds to the evaluation of 10|Jr | differents subsets). Finally, after stopping the meta-heuristic search, the problem (10) is solved by using non-aggregated data for the best subset Jr found, and then the optimal weights are determined. Even though global optimality is not guaranteed using the proposed procedure, its outputs met most of the desired characteristics in a reduced set of scenarios and preserved the statistical properties of the original set of scenarios.

4

Scenario reduction a priori results

This section presents the application of the method described in section 3 to the full set of hydrological scenarios of the Chilean SIC, and discusses how well the moments of the SSoHSc match the moments of the CSoHSc, before the use of the scenarios in a stochastic optimization application (a priori results). The multivariate distribution, as defined in section 2, has 204 dimensions (it consist on the hydroelectric energy inflows to the system, grouped by bus). Subsets from 1 to 14 elements were selected. To evaluate the results of the scenario reduction method over the multivariate distribution it is convenient to define the root mean squared error (RMSE) for the moments as (11), with an equivalent definition for the error on the sub-diagonal matrix of M. The RMSE can be expressed as a percentage of the total moment (i.e. considering all the observations of each component of X as individuals observations of the same process).

  1

RMSE µ(i) = √ µ(i) − m(i) (11) N The preceding definitions allow to show the scenario reduction results in terms of RMSE versus the subset size, as shown in Fig. 2. As expected, the representation improves as the subsets sizes grow. Fig. 2 also shows the results of selecting scenarios by quantiles of the yearly inflow energy (all

9

Historical data: Inflow [m3 /s] detailed by river, month and year

Partial intra-annual aggregation Pre-grouping inflow vectors in n clusters

Partitioning Around Medoids, pam()

Parallelized meta-heuristic search of optimal subset

Differential Evolution, DEoptim()

Estimated valuation of each test subset (Obj. function)

Barzilai-Borwein spectral methods, BBoptim()

stopping criterion

Weights’ re-calculation using non-aggregated data for the “best” subset

BBoptim()

Optimal subset and weights

Figure 1: Flow chart of the scenario selection algorithm.

10

35 25 20 15 10 5 0

(1)

RMSE(µInEnergy ) [%]

30

Quantile Sel. Moment-matching Sel.

2

4

6

8

10

12

14

8

10

12

14

RMSE(µInEnergy ) [%]

(4)

100 60 20 0

RMSE(Msubdiag InEnergy ) [%] 6

6

8

10

12

14

2

4

6

8

10

12

14

200

4

4

150

2

2

100

14

50

12

0

100

10

150

(2)

60 20 0

8

200

6

100

(3)

4

50

RMSE(µInEnergy ) [%]

2

0

RMSE(µInEnergy ) [%]

Subset Size

Figure 2: Moment RMSE of the SSoHSc with respect to the CSoHSc, for the hydroelectric energy inflow to the SIC. equally likely), which is a simpler and more intuitive method (e.g., for three scenarios it would select the first, second, and third quartiles) and frequently used in capacity expansion planning applications. The moment-matching methodology shows to be highly superior, with differences around 3-10% in the error of the first moment and much higher in high-order moments (even close to 80% in the fourth moment).

5

Scenario reduction a posteriori results

This section applies the selected subsets from section 4 to an operational model of the SIC in order to assess how well they would perform in a concrete stochastic optimization application in the electric power industry. This assessment is conducted by measuring the distortion of the SIC operational results using the SSoHSc when compared to results obtained using the CSoHSc (a posteriori results).

11

In operational planning optimization, the objective function basically consists of the expected generation cost. In capacity expansion planning under uncertainty, the objective function can be formulated (in a very simplified way) as in (12), where z correspond to the investment decisions and ω to the stochastic energy inflow. It is clear from the objective functions of both optimization problems that proper representation of the variability in the generation cost under hydro uncertainty is paramount in both cases. Thus, it is very important that the generation cost probability distribution is well represented with the subset of selected hydrological scenarios. min Investment Cost(z) + Eω [Gen Cost(ω, z)] z

(12)

For evaluating the scenario reduction performance in the operation variables of the Chilean SIC, a model of this power system was built based on public information [1] [19], resulting in a model that includes 74 buses and 180 generators. The model database was simulated with PLEXOS [20]. PLEXOS is a market modeling software capable of optimizing unit commitment and electric load dispatch in hydrothermal systems. The resulting mixed-integer linear programs are then solved using XPress [21]. The SIC model was used to simulate one year of operation (2011) for each one of the 48 years of full hydrological data3 (1961-2008) independently. Fig. 3 shows how the monthly generation cost behaves in terms of variations in the energy inflows. Each month of the year is plotted with a different marker, so it is easy to see that under the same conditions the generation costs decreases in a non-linear way with the energy inflows. If we imagine a marginal probability distribution for the energy inflow in Fig. 3, it will be transformed or distorted when passing through the model to obtain the operational cost. In other words, the generation cost could be seen as a transformation of the hydro inflows. Even when the energy inflows are represented properly with the SSoHSc, the power system could distort the representation of the objective function, as observed when comparing Fig. 2 and Fig. 4. However, Fig. 4 shows that despite the distortions in the probability distributions, the moment-matching method allows to represent the expected generation cost with a RMSE (defined similarly to the energy RMSE, but considering only monthly components) below 5% with a SSoHSc of 7 elements or more. Also, in terms of the RMSE of the generation cost, the proposed method behaves better than the quantile selection. 3

For 1961 and 2009 data there are some months missing.

12

200 150 100 50 0

Monthly generation cost [MMUS$]

January February March April May June July August

September October November December 1000

1500

2000

2500

3000

Energy inflow [GWh]

Figure 3: Monthly generation cost as function of the monthly energy inflow in the SIC. At this point, it is important to clarify that the error on the value of the objective function of the operation problem (i.e. the expected generation cost without splitting it into its monthly components) was never superior to 1% in any subset constructed by the moment matching method. Another important output of power system models (especially when dealing with transmission) are flows on the lines and prices on certain nodes. In that sense, Fig. 5 shows that the power flow duration curve of the southnorth transmission corridor (the one that carries most of the yearly hydroelectric energy to the load center of the SIC) can be approximated properly by a SSoHSc of 11 elements, but if the interest is on the peak of the curve it suffices to use SSoHSc of two or five elements (selected with the momentmatching method). The system prices are also well approximated by using the SSoHSc. Fig. 6 shows the price duration curve in the reference node of the SIC. Prices can not be approximated well with a SSoHSc of only two elements, mainly because the two selected hydro scenarios are not capable of stressing the system enough to make the price go above the marginal cost of the natural gas generation, as it would happen if a drier hydro scenario was included in the SSoHSc.

13

35 25 20 15 10 0

5

(1)

RMSE(µOpCost ) [%]

30

Quantile Sel. Moment-matching Sel.

2

4

6

8

10

12

14

8

10

12

14

80 100

(4)

80 60 40 20 0

RMSE(Msubdiag OpCost ) [%] 6

6

8

10

12

14

2

4

6

8

10

12

14

80 100

4

4

60

2

2

40

14

20

12

RMSE(µOpCost ) [%]

10

0

80 60

8

60

(2)

40 20 0

6

40

(3)

4

20

RMSE(µOpCost ) [%]

2

0

RMSE(µOpCost ) [%]

Subset Size

Figure 4: Moment RMSE of the SSoHSc with respect to the CSoHSc, for the generation cost of the SIC.

6

Conclusion

The proposed methodology for selecting a SSoHSc based on the momentmatching method accomplishes the goal of representing the CSoHSc adequately in terms of the multivariate distribution moments and covariances. The selected SSoHSc also showed to be capable of approximating the probability distribution of the objective function of the operational planning problem. Subsets with six or more elements exhibited errors in the expected value under 5% in RMSE, and under 1% in the total value of the generation cost. Other variables of the system, such as prices and flows, are also well represented in terms of their duration curves by using a SSoHSc selected with the proposed method. The results provide us with the necessary confidence to use the selected SSoHSc in further work. Directions of further research include the improvement of this method by using PCD and transformations, so that the selection can be made over 14

1250 1000 750 500 250 0

Full Collection Coll. 2 elems. Coll. 5 elems Coll. 11 elems

-250

“Charrua 500”→“Ancoa 500” Flow [MW]

security limit

0

20

40

60

80

100

Time [%]

Figure 5: Flow duration curve in the corridor “Charrua-Ancoa”, which carries most of the injected hydro energy to the load center of the SIC. Time is in percent of the 8760 hours of one year. auxiliary variables linearly related to the generation cost. Also, the SSoHSc selected with the method are currently being used in different electric power industry applications. Particularly, they are being used for solving the capacity expansion problem in the Chilean SIC using stochastic mixed-integer programming, since using the CSoHSc for this problem is impractical due to the size of the mathematical optimization problem.

References [1] C. N. de Energ´ıa, Fijaci´on de Precios de Nudo Abril de 2011, Sistema Intercontectado Central (SIC): Informe T´ecnico Definitivo, Santiago, Chile (abril 2011). [2] P. W. Glynn, M. Avriel, G. B. Dantzig, Decomposition and parallel processing for large-scale electric power system planning under uncertainty, in: Proceedings of the Workshop on Resource Planning under Uncertainty, Stanford University, 1989, pp. 9–34. [3] I. Aravena, R. C´ ardenas, E. Gil, V. Hinojosa, Capacity expansion planning of the Chilean SIC using stochastic programming, Final Report, report of the DIE UTFSM – Transelec RET I+D Project, Valpara´ıso, Chile (nov 2012). 15

250 200 150 100 50 0

“Quillota 220” Price [US$/MWh]

Full Collection Coll. 2 elems. Coll. 5 elems Coll. 11 elems

0

20

40

60

80

100

Time [%]

Figure 6: Price duration curve for the reference node of the SIC. Time is in percent of the 8760 hours of one year. [4] M. Dyer, L. Stougie, Computational complexity of stochastic programming problems, Mathematical Programming 106 (3) (2006) 423–432. [5] M. H. van der Vlerk, Stochastic mixed-integer programming, in: 12th International Conference on Stochastic Programming (SP XII 2010), Pre-conference Workshop, Halifax, Canada, 2010. [6] T. K. Kristoffersen, Stochastic programming with applications to power systems, Ph.D. thesis, Department of Operations Research, University of Aarhus, Aarhus, Denmark (jul 2007). [7] J. da Costa, G. de Oliveira, L. Legey, Reduced scenario tree generation for mid-term hydrothermal operation planning, in: International Conference on Probabilistic Methods Applied to Power Systems (PMAPS), 2006, pp. 1 –7. doi:10.1109/PMAPS.2006.360199. [8] M. Kaut, S. W. Wallace, Evaluation of scenario-generation methods for stochastic programming, no. 14 in Stochastic Programming E-Print Series, Institut fr Mathematik, 2003. [9] J. Dupaˇcov´ a, N. Gr¨ ove-Kuska, W. R¨omisch, Scenario reduction in stochastic programming: An approach using probability metrics, Mathematical Programming 95 (2003) 493–511. [10] H. Heitsch, W. R¨ omisch, A note on scenario reduction for two-stage stochastic programs, Operations Research Letters 35 (6) (2007) 731– 738.

16

[11] K. Høyland, M. Kaut, S. W. Wallace, A heuristic for moment-matching scenario generation, Computational Optimization and Applications 24 (2-3) (2003) 169–185. [12] Z. Chen, D. Xu, Knowledge-based scenario tree generation methods and application in multiperiod portfolio selection problem, Applied Stochastic Models in Business and Industry (2013) n/a– n/adoi:10.1002/asmb.1970. URL http://dx.doi.org/10.1002/asmb.1970 [13] R Development Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0 (2012). URL http://www.R-project.org/ [14] D. Aloise, P. Hansen, L. Liberti, An improved column generation algorithm for minimum sum-of-squares clustering, Mathematical Programming 131 (1-2) (2012) 195–220. [15] K. Mullen, D. Ardia, D. Gil, D. Windover, J. Cline, DEoptim: An R package for global optimization by differential evolution, Journal of Statistical Software 40 (6) (2011) 1–26. URL http://www.jstatsoft.org/v40/i06/ [16] R. Varadhan, P. Gilbert, BB: An R package for solving a large system of nonlinear equations and for optimizing a high-dimensional nonlinear objective function, Journal of Statistical Software 32 (4) (2009) 1–26. URL http://www.jstatsoft.org/v32/i04/ [17] K. Pearson, On Lines and Planes of Closest Fit to Systems of Points in Space, Philosophical Magazine 2 (11) (1901) 559–572. [18] A. Reynolds, G. Richards, B. de la Iglesia, V. Rayward-Smith, Clustering rules: A comparison of partitioning and hierarchical clustering algorithms, Journal of Mathematical Modelling and Algorithms 5 (2006) 475–504. [19] D. de Peajes CDEC-SIC, Revisi´on 2011 Estudio de Transmisi´on Troncal, Cuatrienio 2011-2014, Santiago, Chile (junio 2011). [20] Energy Exemplar, PLEXOS for Power Systems–Power Market Simulation and Analysis Software [computer software], december 2013. URL http://www.energyexemplar.com/. 17

[21] B. Daniel, Xpress-Optimizer Reference Manual, Fair Isaac Corporation, Leamington Spa, Warwickshire, UK (jun 2009).

18