University of Southern California Dissertations and Theses

The next generation of power-system operations: modeling and optimization innovations to mitigate renewable uncertainty

(USC Thesis Other)

The next generation of power-system operations: modeling and optimization innovations to mitigate renewable uncertainty

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content The Next Generation of Power-System Operations: Modeling and Optimization Innovations to Mitigate Renewable Uncertainty Semih Atakan August, 2018 University of Southern California Epstein Department of Industrial and Systems Engineering Submitted to the faculty of the USC Graduate School in partial fulllment of the requirements for the degree of Doctor of Philosophy in Industrial and Systems Engineering Dissertation Committee: Suvrajeet Sen, Advisor Julia L. Higle Rahul Jain Abstract The momentum of renewable energy is transforming the power grids across the world. A markedly increasing share of electricity-demand is sourced by solar and wind resources, which are volatile and dicult-to-store, therefore causing concern regarding the reliability of power-system operations. This thesis adopts a multi-faceted approach on expanding our knowledge of decision-making tools which could eventually be benecial to mitigate this uncertainty, through improved optimization of day-to-day operations. We begin with the so-called unit-commitment problem and propose a novel formulation that is tailored for problems with longer planning horizons. With its strong polyhedral properties, our formulation can deliver optimal solutions to week-long and industry-strength problems in nearly 10 minutes, using an ordinary desktop computer. Furthermore, our research shows that considering week-long planning horizons, in lieu of the industry-standard day- long ones, captures demand trends better, and promises improved operating-costs. Likewise, this formulation can accommodate the intra-week trends in solar and wind resources, and may produce decisions with better foresight. Our second contribution is a methodological one that aims to address uncertainty in general mixed-integer formulations. In particular, we consider the well-known Progressive Hedging (PH) algorithm from the stochastic convex optimization literature, and extend it for solving multi-stage stochastic mixed-integer convex programs (SMICPs). Our framework for using PH is guaranteed to compute optimal solutions for SMICPs, which contrasts with most previous extensions of PH, as they have been implemented without convergence guarantees. More importantly, we demonstrate the eectiveness of our framework through computational experiments, and observe signicant improvements over the performances of benchmark algorithms. In particular, our algorithms can deliver optimal solutions to certain instances in nearly 30 seconds, whereas a decade ago, multiple hours of computational time had not been enough to identify an optimal decision. Finally, we turn our attention to renewable energy, and assess the risks of increasing levels of renewable- integration on power grids. To this end, we develop an experimental framework that resembles the daily operations of power-management authorities. We also use this framework to assess the use of deterministic and stochastic optimization methodologies for dierent components of power-system planning processes. Our experiments provide quantitative support to the notion that introducing signicant renewable resources into power systems requires a systems-view of supply and demand uncertainty, the transmission capacity of the grid, and just as importantly, the operations planning strategy used. With regards to the latter point, we demonstrate that stochastic programming can be a valuable decision-making tool for maintaining the reliability and cost-eectiveness of the grid. i Acknowledgements My family will almost surely never read the entirety of this dissertation. Yet, I was able to pursue a doctoral degree and produce this dissertation thanks to their moral and nancial support. I am grateful to them for giving me this opportunity. Throughout my studies, I was surrounded by many loving friends who have always managed to keep my spirits high. Among those, I name here only the ones that have understood what my thesis is all about. I am grateful to Miju Ahn and Yifan Liu for their great friendships that made my time at USC enjoyable, and to Harsha Gangammanavar for both his friendship and awesome company in many of our vacations together. Additionally, I am thankful to G ok ce Kahvecio glu for giving me the much-needed emotional support in all times of diculty and for her continuing friendship. Several people have taken the time to help me out during my studies. Among those, I appreciate Dr. Julia L. Higle, Dr. Rahul Jain, and my advisor for being a part of my dissertation and qualifying exam committees, and Dr. Maged Dessouky, and Dr. Meisam Razaviyayn for being a part of my qualifying exam committee. In particular, Dr. Higle's feedbacks have played an important role in clarifying the exposition of this dissertation. Furthermore, the ISE department faculty and sta have always been helpful and supporting, and I am grateful for it. Certain chapters of this thesis rely on the research conducted in collaboration with Dr. Guglielmo Lulli and Dr. Harsha Gangammanavar, along with my advisor. I would like to thank them for their contributions and guidance that played a signicant role in forming my dissertation. Additionally, a version of Chapter 3 will appear in a journal's special issue honoring Dr. Maarten van der Vlerk. It is my pleasure to honor his life-long contributions to the stochastic programming literature through this dissertation as well. My master advisors, whom I thought I had left in Turkey, have, in fact, always been with me throughout my studies. I am grateful to Dr. Kerem B ulb ul and Dr. Nilay Noyan for their ceaseless support and friendship, and their huge impact on the course of my career. Finally, a few words must be spared for the person who made this dissertation happen. I am grateful to my advisor, Dr. Suvrajeet Sen, for his guidance and support that played a tremendous role in shaping this dissertation. Many of the ideas presented in here bear his ngerprints, which I had tried to perfect in my growing capacity. It was his hard-earned grants that allowed me to participate in over 10 professional meetings that deepened my academic abilities, and I appreciate this. More importantly, I am indebted to Dr. Sen for his never-ending personal support and encouragements, together with his friendly and patient attitude. His impact on my personality and career will remain to be signicant throughout my life. This thesis, truly, would not have been possible without him. Support for this research comes through the grants by the National Science Foundation (ECCS-1548846) and the Air Force Oce of Scientic Research (FA9550-13-1-0015). ii Contents Abstract i Acknowledgements ii 1 Introduction 1 1.1 Power-System Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 The Unit Commitment Problem and Extensions . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Background in Mathematical Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Mixed-Integer Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 Stochastic Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 A State-Transition MIP Formulation for the UC Problem 8 2.1 Formulating Production and Operations in Unit Commitment . . . . . . . . . . . . . . . . . . . . . . 8 2.1.1 Formulating the Components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.2 Benchmark Formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 The State-Transition Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Insights on the Polyhedral Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3 Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3 The PH-BAB Framework 26 3.1 Progressive Hedging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 General Framework of PH with Branch and Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3 PH-BAB for a Special Class of Two-Stage Stochastic Programs . . . . . . . . . . . . . . . . . . . . . 37 3.4 Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.1 Experimental Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.2 Performance Comparisons on Linear SMIP Instances . . . . . . . . . . . . . . . . . . . . . . . 43 3.4.3 Sensitivity of PH-based Algorithms to . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.4.4 Performance Comparisons on Quadratic SMIP Instances . . . . . . . . . . . . . . . . . . . . . 45 3.4.5 Parallel Computations in the PH-BAB Framework . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.6 Performance of PH-BAB Beyond Specially-Structured Problems . . . . . . . . . . . . . . . . 48 3.5 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4 Taming the Duck: Can Stochastic Optimization Help? 51 4.1 The Quandary of Renewable Energy Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Models for Power-System Operations Planning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.1 Modeling Daily Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 iii CONTENTS iv 4.2.2 Optimizing Commitment and Generation Schedules . . . . . . . . . . . . . . . . . . . . . . . 55 4.2.3 Predicting Solar and Wind Outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.4 Simulating Power-System Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3 Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.4 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5 Conclusion 67 A Data Generation for Deterministic UC Problems 69 A.1 Synthetic Instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 A.2 Realistic Instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 B Stochastic Server Location Problem SMIP Formulation 71 Acronyms 73 Bibliography 74 Chapter 1 Introduction 1.1. Power-System Operations Over the last century, electricity has become the most widely consumed form of energy source in the world. Modern societies maintain ceaseless access to electricity, where the consumers are reminded of its presence, only, thanks to their bills. Similar to its spread, the eciency and reliability of electricity production and distribution have made remarkable progress. To give an example, in the 40 years following 1929, the price of electricity has plummeted from 60¢ to 10¢ per kilowatt-hours (in 2005 dollars) for residential customers in the U.S., and the gure remains approximately the same today (FERC, 2015). Despite this progress, power-system operations are yet to be deemed as perfect. The massive size of the industry often hinders the use of state-of-the-art decision technologies, as these technologies, in general, may oer limited scalability. Moreover, environmental trends, evolving customer needs, and generator innovations have been transforming the grid, and leading to new decision-making challenges. It is good fortune that such challenges can usually be surmounted with improved decision-models and algorithmic novelties, which is what, in its capacity, this thesis intends. We begin with a general overview of operations in modern power systems. Much like any other physical product, electricity can be quantied, priced, and traded among multiple parties. Most developed countries have liberal electricity markets where generator companies present their bids for producing electricity, and utilities can buy power in real-time to cover the demand in their jurisdictions. These markets are overseen by entities, such as the independent system operators (ISOs), which ensure that the supply of electricity meets the demand at market- clearing prices. Consequently, in industries managed by ISOs or similar authorities, the price of electricity, the demand, and the supply, are all driven by microeconomic principles. Unlike most physical products, the demand for electricity is, in general, signicantly inelastic. A typical con- sumer does not observe the real-time price of electricity, and even if she did, may not be quick to or willing to respond. Technological innovations for designing smarter grids | whose components can respond to electricity prices instantaneously | are still in their infancy. Similarly, the storage of electricity, which can also mitigate demand inelasticity, is only achievable at small scales. Therefore, electricity, in most practical situations, must be produced and consumed contemporaneously by the respective participants of the market. This suggests that if the supply and demand are to match precisely, it is the responsibility of the suppliers (i.e., the generators) to chase the demand, in order to prevent blackouts. However, the demand for electricity is neither deterministic nor stationary. A typical load pattern exhibits high variations with respect to time of the day, day of the week, season, and most importantly, the weather. Therefore, accommodating this variability requires signicant exibility from the generating units. On the other hand, not all generators can comport with this. Most generators are limited by their physical capabilities which may restrict sudden adjustments in their energy outputs, or force them to remain operational, needlessly long, after being put into service | a nuclear reactor cannot be expected to consume uranium intermittently. Other load-following 1 CHAPTER 1. INTRODUCTION 2 generators are much more exible, but their production costs can be higher (e.g., natural gas turbines), or their output may depend on the availability of natural resources (e.g., hydroelectric power plants). Given the scale of most power systems, the diversity in the physical attributes of generators leads to massive trade-os, which are extremely dicult to assess, without the use sophisticated decision mechanisms. Most modern power industries analyze and improve their operations through the use of mathematical models. For instance, these industries routinely solve the unit commitment and the economic dispatch problems, where the former problem determines hourly production schedules of generators for the day-ahead, whereas the latter obtains the best production amounts that can guide real-time operations. The mathematical models of these problems must be designed to re ect the actual system as accurately as possible, and at the same time, must be amenable to optimization with reasonable computational eort. These two objectives often con ict, leading to sacrices in model delity for faster decisions, perhaps with relatively poor quality. The interest of the academic community | and this thesis | to the aforementioned problems stems from such trade-os, as they can often times be eliminated with the development of better technology. Contemporary Challenges Over the last decade, governments have set bold targets to incorporate renewable generators | power plants that are fueled by renewable resources like wind or solar | into their portfolio of energy sources. In 2009, the European Commission set a target of 20% renewable integration by the year 2020, with national targets ranging from a low of 10% (in Malta) to a high of 49% (in Sweden). In 2013, Barack H. Obama, as the president of the United States, demanded the same 20% target from all U.S. agencies, and most recently, China has announced its intentions to spend $360 billion through 2020 on renewable energy sources (Forsythe, 2017). Despite the political turmoils that may lie ahead, the growing trend for renewable integration into the grid appears to be permanent (Obama, 2017). Indeed, public support for expanding wind and solar energy is very strong among Americans (Funk and Kennedy, 2016), and overwhelming { exceeding 90% support { across Europeans (Eurobarometer, 2015). Moreover, several ISOs have already restructured their energy markets and operations in order to better accommodate the \erratic" nature of renewables | ultimately, a wind turbine can only produce when the wind blows, and a solar generator cannot operate after sunset. While the environmental benets of renewable energy are not up to debate, the widespread use of such energy resources poses major problems for power-system operations. A popular example that demonstrates the side-eects of renewable integration is the so-called duck curve (or duck-chart), which was coined by California ISO in 2013, and is now routinely picked up by media outlets (see, for instance, The Economist, 2017). This curve depicts the aggregated daily net-load (load minus solar and wind supply) of the grid, facing the prospect of high-levels of solar-based energy, which leads to a growing abundance of supply during daytime. As a consequence, the net load plunges when the sun is up, and surges back during sunset, in a fashion that resembles the \neck" of a duck, thereby earning the moniker (see California ISO (2016) for a detailed description; see, also, Denholm et al. (2008) for proto-depictions of the duck curve). This phenomenon already upsets the operations by causing costly ramps in production by fossil-fueled generators, occasionally leading to \negative" prices in the electricity market (see Denholm et al., 2015, and the references therein). If not addressed, the situation will surely worsen with further solar penetration. Notable investments have been made to ameliorate the negative eects of renewable energy on power-system operations. For instance, in 2008, the New York ISO established a centralized wind-forecasting system to improve the forecasting of wind availability and timing of wind-powered generation (New York ISO, 2016). Such investments are clearly essential for maintaining a healthy and cost-ecient grid, but by themselves, they may not be able to achieve such objectives. Indeed, novel decision models must be adopted, which can harness the information supplied by these systems, and produce decisions that are less susceptible to exogenous factors (e.g., the weather) aecting the industry. CHAPTER 1. INTRODUCTION 3 1.1.1. The Unit Commitment Problem and Extensions The unit commitment (UC) problem has a particular importance to this thesis, and therefore, deserves a detailed treatment. Every day, regional electricity networks deliver millions of kilowatt-hours of energy from generating units to consumers, and these production requirements signicantly vary with respect to the day of the week, and hour. As a result, ecient scheduling of electricity production continues to attract signicant attention from both industry and academia in the form of the so-called UC Problem. Such a model must recommend which generators to use, and how much they should produce so that demand over an hourly-discretized planning horizon is met, while obeying certain operating rules, maintenance schedules, and in some instances, transmission capacity requirements. The goal is to obtain the most cost-eective operating schedule over a large set of generators, and ensure that the demand can be completely fullled during the planning horizon. With the liberalization of the energy industry, the introduction of energy market concepts, and privatization, the role of the UC models have changed (Hobbs et al., 2001). UC problems are solved by a variety of constituents of the electric-power industry, each for a dierent purpose. A specic utility may solve the UC problem for the purpose of planning production within their delivery area. Their costs also help them place bids in the ISO exchange market. On the other hand, the ISO solves its UC model to decide which bids to accept, and to determine prices that will be paid to the suppliers. These problems are much larger because they involve multiple suppliers, and with a number of generating units in the order of thousands. Providing high-quality solutions to these realistic problems is computationally demanding, but has the potential for signicant reduction in costs. A 2011 report of the Federal Energy Regulatory Commission (FERC) suggests that savings approaching $100 million per annum can be expected by replacing heuristics with methods that seek optimal solutions (see O'Neill et al. (2011)). Consequently, developing solution methods that can achieve high-quality solutions in short amounts of time has been the focus of signicant research over the last several decades. Many optimization methods have been proposed to solve the UC problem. For example, we mention branch- and-bound methods (Turgeon, 1977), dynamic programming approaches (Pang et al., 1981), Lagrangian-relaxation methods (Muckstadt and Koenig, 1977; Bard, 1988), and unit decommitment (Tseng et al., 2000), among others. For a detailed review, the reader is referred to (Padhy, 2004; Saravanan et al., 2013). In recent years, mixed- integer programming (MIP) has emerged as a popular tool for solving UC problems. A discussion on its merits and drawbacks (with respect to Lagrangian relaxation) were given in Li and Shahidehpour (2005). The popularity of MIP led to signicant advances in the mathematical description of specic features of the UC problem. For instance, Rajan and Takriti (2005) have identied the convex hull for a minimum up/downtime polytope, Gentile et al. (2016) provided the convex hull of generation limits and this minimum up/downtime polytope, Ostrowski et al. (2012) provided strengthened inequalities of a ramping polytope, and Damc-Kurt et al. (2015) provided the convex hull of a two-period ramping polytope. Incorporating Uncertainty Due to the complexity of UC problems, the traditional approach for addressing uncertainty has been, in essence, over-producing, which creates a reserve of electricity that is to be used as a buer against unexpected volatility. Over the last decades, developments in stochastic modeling and optimization have allowed the uncertainty to be explicitly handled within the decision process. While the reserve requirements are still a part of power-system operations, better approaches, such as stochastic programming, have allowed smaller reserve amounts, thereby reducing the operational costs, and at the same time, maintaining reliability (see, for instance, Ruiz et al., 2009). Early applications of stochastic optimization on UC problems mainly focused on the uncertainty in demand (or price), rather than supply, as renewable energy sources had not yet become as prevalent as they are today. In order to address demand uncertainty, Carpentier et al. (1996) and Takriti et al. (1996) consider multi-stage stochastic mixed-integer linear programs (SMILPs), relying on a decomposition algorithm of Cohen and Zhu (1984), and the Progressive Hedging algorithm of Rockafellar and Wets (1991), respectively. Shiina and Birge (2004) attempt to solve a similar multi-stage SMILP with a column-generation algorithm. Takriti and Birge (2000) incorporate power trading and price uncertainty into their modeling, and solve the resulting problem through a combination of CHAPTER 1. INTRODUCTION 4 Lagrangian-relaxation and Benders'-decomposition methods. Finally, a chance-constrained optimization approach is given by Ozturk et al. (2004). With the spread of renewable energy, supply uncertainty has become the primary concern of the industry, and thereafter, received signicant attention from the academic community. For mitigating wind uncertainty, Wang et al. (2008) present a security-constrained UC problem, and Papavasiliou et al. (2011) investigate improved reserve- commitment through a two-stage SMILP, solved with the dual decomposition algorithm of Care and Schultz (1999). Bertsimas et al. (2013) adopt a robust optimization approach that can be solved with a Benders'-decomposition scheme. Through the use of stochastic UC models, the impact of wind volatility on electricity markets has been investigated by Barth et al. (2006); Abrell and Kunz (2015); Zavala et al. (2017). In addition to supply uncertainty, other sources of stochasticity | such as generator and transmission failures, and/or demand uncertainty | have also been considered in Papavasiliou and Oren (2013); Zheng et al. (2013); Cheung et al. (2015). 1.2. Background in Mathematical Optimization The decision models and solution algorithms that are considered in this thesis will be based on two branches of mathematical optimization, namely mixed-integer programming and stochastic programming. In what follows, we provide preliminary background on both of these subelds. For more detailed treatments, we refer the readers to Nemhauser and Wolsey (1999) for mixed-integer programming, and Kall and Wallace (1994); Birge and Louveaux (1997) for stochastic programming. 1.2.1. Mixed-Integer Programming A wide variety of operational and nancial problems are modeled as mixed-integer programs (MIPs) where a certain subset of the decision variables are restricted to assume integer values. Examples include scheduling, production and inventory planning, transportation, portfolio optimization, among many others. The most typical forms of MIPs assume linear objectives and polyhedral feasible sets, thereby labeled as mixed-integer linear programs (MILPs). A mathematical statement of these problems can be given as follows: min x c | xjAx =b; x 0; x2R nr Z nz : (MILP) The above problem includes n variables, n z of which are restricted to integer values, and the remaining n r are continuous. These problems are inherently dicult to handle as the decisions must be determined over a discrete (therefore, non-convex) feasible set. A broader class of MIPs is obtained when linearity of the objective function and constraints is replaced with convexity of these objects. This leads to mixed-integer convex programs (MICPs) which can be stated as min x f(x)jx2C; x2R nr Z nz : (MICP) Above, we assume a closed and convex set C, and a (closed) convex function f. Solving MICPs is naturally more challenging than solving MILPs, given that ordinary convex programs, in general, are relatively more dicult to solve in comparison to linear programs (LPs). However, eective algorithms for special cases (e.g., mixed-integer quadratic programs) are known to exist and used for a wide variety of purposes. The state-of-the-art solvers that address MIPs typically rely on branch-and-bound strategies. Over the last decades, the inclusion of cutting planes and heuristics into these solvers has improved their performances remarkably (see Bixby, 2002). However, it is well known that MIPs are, in general,NP-hard. Despite all the advances, solving problems of realistic sizes | involving millions of variables and constraints | may still not be trivial. As a result, there has been signicant attention on developing well-crafted MIP formulations, which can induce high quality solutions, in short amounts of time. For instance, a natural strategy is to formulate the constraints of the problem, CHAPTER 1. INTRODUCTION 5 i.e.,C, in such a way that it resembles | as closely as possible | the convex hull of the feasible region of the problem, denoted as conv(C\ (R nr Z nz )). On the other hand, a formulation should not be ooded with constraints, as it would reduce the pace at which the relaxations in a branch-and-bound tree can be processed. Maintaining a delicate balance between these two ends is what makes an MIP formulation amenable to producing fast and high-quality solutions. 1.2.2. Stochastic Programming An important but challenging venue for researchers and practitioners is the incorporation of data uncertainty into the mathematical formulations of optimization problems. In real life, much of the input data is usually not known with certainty at the time the decisions are made. Ignoring the inherent uncertainty in the systems (for instance, by using the expected values of the random parameters) can lead to unfavorable decisions, which, in general, underperform in real-life situations. A general methodology to address uncertainty is to consider a sample of scenarios as a proxy for the random parameters, and optimize the decisions based on these outcomes. This leads to stochastic programming formulations, which are notorious even without the presence of integrality requirements on the decisions. The reason is, the large numbers of samples that are necessary for statistical accuracy, often times lead to massive scale problems, which, naturally, deteriorate the computational performance of the underlying solvers. Preliminaries In a stochastic program, the considered scenarios can either be formed by extrapolating from the past data, or more rigorously, by assuming the existence of random variables that govern the modeled system. We represent the uncertainty with the random vector , whose realizations s ; s = 1:::S, re ect the state of the problem under scenarios. We denote the set of scenariosf1:::Sg withS. The random variables (i.e., components of) are dened over the probability space ( ;F;P). The probability of a scenario is then denoted with p s =P ((!) = s ) where p s > 0;8s, and P S s=1 p s = 1. A nite-dimensional problem is posed by assuming either the sample space has nite support, or letting S represent the nite number of sampled scenarios from an arbitrary distribution. The goal of a stochastic program is to nd the best policy for the underlying stochastic optimization problem. A policy X() should be interpreted as a random decision-vector whose response to scenario s is a decision vector X( s ) = x s . We nd it convenient to represent a policy with the shorthand mapping X : s! x s (more formally, the mapping X should be dened as X :S!R n , but we adopt the former to make the outcomes of the mappings more explicit). A policy X can be partitioned into T components (X 1 ;:::;X T ) where X t should be interpreted as the segment of the policy to be followed in time period t. Likewise, ltrationF of the probability space is partitioned into periods asF 1 :::F T , whereF t F t 0 for t 0 t. A logical requirement is that an implementable policy in an arbitrary stage should have no foresight on the exact realizations of future events. Consequently, for any set of scenarios FF t which are observationally indistinguishable at time t, we must have x st ;8s2F , equal to each other. For instance, at the initial time period, a policy-maker cannot be anticipative of the future, and therefore must take a single action that is not contingent on any specic scenario (i.e., x s1 = x;8s = 1:::S). In the stochastic programming literature, such restrictions are commonly referred to as the nonanticipativity constraints. These constraints will be collectively represented by the set N . Accordingly, we can ensure that a policy is nonanticipative with the restriction X2N . Stochastic Convex Programming We begin with the mathematical statement of stochastic programs which involves no integrality restrictions on its feasible set. A broad and relatively well-understood family of problems that suits this assumption is the stochastic convex program (SCP), which can be dened as follows: min X E[F (X) ] subject to X2N\ S s=1 C s : (SCP) CHAPTER 1. INTRODUCTION 6 Above, the mapping F (X) : s! f s (x s ) embodies the scenario objectives functions f s (x s ). The expectation of F (X) is dened asE[F (X) ] = P S s=1 p s f s (x s ). For each subproblem, we haveC s as a closed and convex set, andf s as a (closed) convex objective function. The nonanticipativity setN admits a linear description, therefore the SCP also encapsulates its special case, the stochastic linear program (SLP), where f s are linear and C s are polyhedral, for every scenario s. Observe that the restriction X2N prevents the problem from decomposing into individual scenario subproblems, immediately. Eective solution methodologies for stochastic programming rely heavily on decomposition ideas, where the decomposition may either be performed with respect to the scenarios, or the time periods (a.k.a., stages) in the decision model. The L-shaped method (Van Slyke and Wets, 1969) for two-stage SLPs, and the nested Benders' decomposition (Louveaux, 1980; Birge, 1985) for multi-stage SLPs can be considered as representatives of the rst category. Examples within the context of scenario decomposition involve the Progressive Hedging algorithm of Rockafellar and Wets (1991) for multi-stage SCPs, and the diagonal quadratic approximation of Mulvey and Ruszczy nski (1995) for multi-stage SLPs. Stochastic Mixed-Integer Convex Programming When a subset of decisions in an SCP is required to satisfy integrality restrictions, the resulting problem becomes a stochastic mixed-integer convex program (SMICP), and can be stated as follows: min X E[F (X) ] subject to X2N\ S s=1 C s ; x s 2R nr Z nz ;8s = 1:::S: (SMICP) For every scenario subproblem s, the decision vector x s has n components, n z of which are restricted to an integer lattice, and the remaining n r are continuous. The rest of the denitions are analogous to the SCPs. It is worthy to mention that when f s (x s ) is linear and minimized over a polyhedral set (which would lead to an MILP), the SMICP reduces to an SMILP. For SMICPs, convergence to a global optimal solution in a computationally-eective manner is hard to achieve, not only due to high numbers of coupled subproblems, but also due to the non-convex nature of the feasible set. In most cases, global optimality is ensured | under the assumption of linearity | via branch-and-bound- or cutting-plane-based procedures. For instance, for multi-stage SMILPs, the dual decomposition of Care and Schultz (1999) and the branch-and-price algorithm of Lulli and Sen (2004) use scenario decomposition within a branch-and-bound framework. For two-stage SMILPs, cutting-planes can be used for the proper convexication of the second-stage value functions. For instance, the algorithms of Qi and Sen (2016) and Ralphs and Hassanzadeh (2014) are applicable to a wide variety of problems. Further assumptions (such as, pure or mixed integers, binary or general integers) on the individual stages lead to computationally faster procedures (see Laporte and Louveaux, 1993; Care and Tind, 1998; Sen and Higle, 2005; Sen and Sherali, 2006; Gade et al., 2014; Zhang and K u c ukyavuz, 2014; Angulo et al., 2016). A very informative table that summarizes the characteristics of several algorithms can be found in Ralphs and Hassanzadeh (2014) (see also Trapp et al. (2013)). 1.3. Contributions The main purpose of this thesis is to put together a promising set of tools that address major decision-making challenges faced by the power industry. In light of the transformations that power systems are going through, this thesis adopted the theme of mitigating uncertainty in the grids of the future. To this end, we oer both applied and methodological contributions that are tailored, specically, to the domains of power-system operations planning and stochastic mixed-integer programming. More specic contributions are summarized as follows: In Chapter 2, we analyze the well-regarded MIP formulations of a deterministic UC problem, and propose a novel formulation that can tame its computational complexity. The new formulation has many desirable properties, such as naturally-appearing facet-dening constraints that lessen the need for additional valid inequalities. Furthermore, CHAPTER 1. INTRODUCTION 7 the new formulations demonstrates superior performance for problem instances with longer planning horizons. This chapter also demonstrates the added value of considering longer planning horizons in UC problems, and suggests our formulation as an eective tool for reaping the exemplied prospects. In Chapter 3, we turn our attention to the theory of stochastic programming, and develop a new framework for optimizing multi-stage SMICPs. The suggested framework is a decomposition scheme capable of providing optimality guarantees on the produced decisions. Unlike most existing algorithms, the theoretical foundation of our framework supports the use of arbitrary convex objectives and feasible sets, therefore, greatly expands its potential of applicability to a wide range of domains. We demonstrate the eciency of the proposed framework through well-known benchmarks, and their extensions. Our experiments also reveal the decade-long progress made in (stochastic) mixed-integer programming: Problems which could not be solved to optimality in two hours (Ntaimo and Sen, 2005) can now be optimized in half a minute. In summary, these two chapters aim to improve our understanding of deterministic UC formulations and SMICP methodology, respectively. Specic remarks at the end of each chapter suggest how these developments can serve as tools to mitigate renewable uncertainty. The next chapter of this thesis will have a slightly dierent avor and a unifying intention: In Chapter 4, we make a direct attempt to address the issues concerning renewable integration into power grids. In particular, we focus on the uncertainty associated with the renewable energy resources, such as solar and wind. We begin with setting up an experimental framework for modeling major planning processes in (daily) power- system operations. Using this framework, we analyze the impact of renewable integration in terms of operating costs, reliability, and environmental metrics. Furthermore, our study evaluates the value of stochastic programming on common planning problems, and make comparisons to the deterministic approaches used in the industry. Finally, we conclude this dissertation in Chapter 5, with a summary of our ndings, and suggestions for future research directions. Chapter 2 A State-Transition MIP Formulation for the UC Problem In this chapter, we study the MIP formulations of the UC problem. Our focus will be on the components of the problem that are of critical interest when forming a realistic model. More advanced issues, such as the optimization of the transmission network and the handling of supply uncertainty, are left out of consideration. These aspects will be addressed in Chapter 4. The contribution of this study is threefold. First, we develop a new formulation using a novel set of state- transition variables. These variables are named after their role in capturing the transition of generator states between consecutive time periods. They introduce a network sub-structure into the UC problem formulation. As observed originally by Pritsker et al. (1969), and more recently by Bertsimas et al. (2014), this structure is much more amenable to LP-based solution methods, leading to faster solution times, without modifying the underlying optimization methodologies. Second, we show that the use of state-transition variables naturally lead to certain facet-dening constraints in our problem-of-interest. This contrasts with the majority of other formulations available in the literature, which often require the addition of further inequalities in order to generate the same strong (facet- dening) valid inequalities. Finally, we perform a computational study to demonstrate the behaviors of the new, and two relatively recent and well regarded benchmark formulations. Our experiments are performed on two classes of instances: (i) synthetic instances used by other academic researchers, and (ii) realistic instances made available by the FERC for computational testing. This study sheds further light on the classes of instances for which the new formulation outperforms the contenders. The organization of this chapter is as follows. In §2.1, we begin with describing typical formulation strategies for our UC problem. In §2.2, we present our formulation, and in §2.2.1, we demonstrate its polyhedral properties. Our computational experiments are presented in §2.3, and details on the experimented instances are given in Appendix A. 2.1. Formulating Production and Operations in Unit Commitment Beginning with the rst MIP formulation of Garver (1962), a plethora of mathematical formulations have been proposed in the literature to solve instances of the UC problem. These have extended the work of Garver in several directions. The mathematical modeling has been enriched with many additional aspects of the problem, and made it more realistic. For instance, several operational and technological restrictions have been included in The research in this chapter was conducted in collaboration with Guglielmo Lulli, PhD, from the Department of Management Science, Lancaster University, Lancaster, United Kingdom (e-mail: g.lulli@lancaster.ac.uk). A version of this chapter has appeared in Atakan et al. (2018). 8 CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 9 the UC formulations, and signicant eorts have been made in improving the formulation of operating costs. Other objectives have also been considered, such as the minimization of no load or turn o costs, maximization of social welfare, or maximization of the prot of generator companies (see Shahidehpour et al. (2002)). In addition to rening the mathematical modeling, a lot of academic research in this eld has been devoted to developing strong MIP formulations, so that real-scale instances of the problem can be handled with advanced MIP solvers. As considered in Ostrowski et al. (2012) and Morales-Espa~ na et al. (2013), the core requirements of a UC problem involves the following components: • minimum and maximum production restrictions, • start up / shut down limits, and ramping restrictions, • minimum up/downtime requirements, • demand and reserve requirements (the latter is modeled as stand-by capacities) • start up costs (dened as step-functions of the generators' idle times) and production costs (dened as piecewise-linear convex functions of the production amounts) In this chapter, we adopt this core setup although other advanced requirements (e.g., transmission, power ow, line ow, and voltage limits, etc) are typically accommodated in realistic applications Baldick (1995). The purpose of our study will be to compare three alternative MIP formulations, which consider the above listed modeling considerations of a UC problem. The notation for these formulations will be introduced gradually throughout this section. A nomenclature is also provided at the end of this section in Table 2.1. To begin our discussion, we present a prototypical MIP formulation of this UC model. Given a set of generators G, and hourly discretized time periodsT , we introduce the following decision variables that are ubiquitous in the literature: x g;t : State variable (1 if g is operational at time t, 0 otherwise), s g;t =z g;t : Start up / shut down variable (1 if g is turned on / o at time t, 0 otherwise), p g;t : Amount of production by g at time t. The state variables are fundamental for the scheduling of generators. The start up / shut down variables are used to formulate the operating costs of the generators, whereas the production variables determine the dispatch amounts. In what follows, a vector of variables of the same type (say x g;t ;8g2G;t2T ) will be typed in bold (say x). The objective of the mathematical formulation is to compute a power generation plan which satises production requirements along with operational constraints. The total operational costs include both the production and the start up costs. Both of these costs are nonlinear, in general. In our prototypical formulation, these costs are represented with the following (nonlinear) functions: F g (), for the start up costs, and V g (), for the production costs. The argument of the cost functions are respectively i- and j-dimensional vectors x g;[t] = (x g;ti :::x g;t ) and p g;[t] = (p g;tj :::p g;t ), for some i 0 and j 0. This indicates the dependency of the costs, respectively, on the past states and the production levels of the generator. A prototypical formulation is given below: min X g2G X t2T F g (x g;[t] ) +V g p g;[t] subject to x g;t x g;t1 =s g;t z g;t 8g2G; t2T; (2.1a) p g;t C g x g;t 8g2G; t2T; (2.1b) p g;t C g x g;t 8g2G; t2T; (2.1c) P g2G p g;t d t 8t2T; (2.1d) (x; s; z; p)2D (2.1e) (x; s; z)2f0; 1g 3jGjjTj ; p 0: (2.1f) Above,d t is the demand for electricity at period t, C g is the maximum generation capacity and C g is the minimum required production amount when the generator is operational. CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 10 The constraints (2.1a) links thes g;t andz g;t variables to the state variables. Indeed, these variables are completely determined once the values of x are known. They are introduced exclusively to capture the transition of a generator between idle and operational states, and to formulate the operating costs of a generator. Constraints (2.1b) and (2.1c) model the lower and upper bounds, respectively. Constraints (2.1d) impose that the production levels must meet the demand for energy at any period. Finally, constraints (2.1e) state that a feasible solution must belong to the polyhedronD, which will be described in detail in the following paragraphs. The integrality requirements on the (s; z) variables can be relaxed without invalidating the formulation. This observation was exploited by Carri on et al. (2006) to formulate the UC problem with integrality restrictions on x alone. However, the assumed benet of using considerably smaller number of binary variables does not necessarily lead to superior computational performance, as observed by Ostrowski et al. (2012). This is due to improved formulations and more robust MIP solvers. Before proceeding further, we dene the following notation which will be used throughout this paper: R g = R g : Ramp-up / ramp-down limit, S g = S g : Start up / shut down limit, UT g =DT g : Minimum uptime / downtime limit. The data of a typical UC instance usually obey the following relations: C g R g S g C g > 0 and C g R g S g C g > 0; 8g2G. We assume UT g 1 and DT g 1 to avoid cases where the generator is simultaneously turned on and o. 2.1.1. Formulating the Components The majority of the UC formulations in the literature can be perceived as extensions of the prototypical formu- lation (2.1). These formulations dier in the way they dene operational constraints, nonlinear objective functions, and production quantities. Nonetheless, the state of a generator (i.e., generator's commitment status and history) is almost always determined by (x; s; z) variables, as rst dened in the seminal work of Garver (1962). In this section, we will identify all the components that go into modeling dierent considerations in a UC problem. In the next section, we will integrate them together into specic formulations, which correspond to the studies of Ostrowski et al. (2012) and Morales-Espa~ na et al. (2013). In what follows, we will allow nonpositive indices of the decision variables to make it explicit that a solution of the problem might depend on the past states of the generators. In the actual implementation, we x such variables to their realized values according to the available past data. Formulating production quantities An integral part of a UC formulation is the description of feasible produc- tion schedules that consider demand and operating reserve requirements along with minimum/maximum production restrictions. In the UC literature, the predominant choice for formulating production amounts has been through the p g;t variables. These variables are intuitive and simple, however, they may inadvertently introduce unneeded complexity into an MIP formulation. In particular, the presence of minimum production restrictions (2.1b) implic- itly alters the continuous nature of these variables, leading to semi-continuous variables that must either equal 0 or lie within the range [ C g ; C g ]. To avoid this, the following variables can be used in lieu of p g;t : p 0 g;t : The production amount beyond C g provided by generator g at time t. This new variable only accounts for the variable portion of production. The xed portion of production (i.e., the minimum production amount) is associated with the state variables as C g x g;t . Combining these, we can form the relationp g;t = (p 0 g;t + C g x g;t ) which can be used to infer the total production amount. Moreover, after applying this substitution throughout the formulation, one can observe that the minimum production constraints are immediately satised. Consequently,jGjjTj constraints in (2.1b) can be omitted from the formulation without sacricing its delity. The idea to treat the minimum and the variable production amounts separately was originally considered in Garver (1962). Over the decades, this idea had been set aside, until the study of Morales-Espa~ na et al. (2013), who CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 11 provided a modern look at this representation. Along with the use of p 0 g;t variables, their study provides tighter descriptions of the production capacity restrictions which will be presented later in this section. We continue with formulating the operating reserves. These requirements are stand-by capacities that must be kept ready to provide for unplanned outages of generating units. System operators use operating reserves to maintain system reliability and to ensure that the supply-demand balance is achieved seamlessly. Baldick (1995) and Arroyo and Conejo (2000) provide examples on how to formulate these requirements. To demonstrate, we introduce the following alternate sets of variables: p g;t : The maximum generation amount that g can supply at time t, r g;t : The generation amount that g can supply at time t for reserve requirements. Both of these variables, by themselves, capture all the necessary information to formulate the reserve requirements. Indeed, the relation between these variables can be described with the equation p g;t = (p g;t +r g;t ). In view of the discussion in this section, we provide two alternatives for formulating production, which re- spectively appeared in the state-of-the-art formulations of Ostrowski et al. (2012) (see constraints (2.2)-(2.3)) and Morales-Espa~ na et al. (2013) (see (2.4)-(2.6)). We begin with the former, where the (p g;t ; p g;t ) variables are used and the following production limits are imposed for all g2G; t2T : C g x g;t p g;t p g;t C g x g;t + ( S g C g )z g;t+1 : (2.2) The above constraints ensure that (i) if the generator is idle (i.e., x g;t = 0), the production amounts cannot be positive, (ii) if the generator is operational (i.e., x g;t = 1), the minimum production requirements and capacity restrictions must be obeyed, and (iii) if the generator is scheduled to be turned o in the next period (i.e.,z g;t+1 = 1), the maximum generation amount cannot exceed the shut down limit. For t =jTj, the z g;t+1 variable on the right- most inequality in (2.2) is assumed to be 0. The p g;t variables must also be smaller than the start up limit S g whenever the generator is turned on at time t, however, this requirement will already be satised by the ramping inequalities presented in the next section. Finally, letting t be the required reserve amount at timet, the following constraints make sure that the reserve requirements are fullled: X g2G p g;t d t + t 8t2T: (2.3) As an alternative for the above formulation of production, Morales-Espa~ na et al. (2013) utilize the (p 0 g;t ; r g;t ) variables. As previously stated, the minimum production constraints are redundant when p 0 g;t variables are in use, therefore omitted. The following constraints provide a tight description of production limits: p 0 g;t +r g;t ( C g C g )x g;t ( C g S g )s g;t max( S g S g ; 0)z g;t+1 ; 8g2G; t2T; where UT g = 1; (2.4a) p 0 g;t +r g;t ( C g C g )x g;t ( C g S g )z g;t+1 max( S g S g ; 0)s g;t ; 8g2G; t2T; where UT g = 1; (2.4b) p 0 g;t +r g;t ( C g C g )x g;t ( C g S g )s g;t ( C g S g )z g;t+1 ; 8g2G; t2T; where UT g > 1: (2.4c) In the above inequalities, z g;t+1 = 0 whenever t =jTj. As shown in Gentile et al. (2016), these inequalities are facets of the convex hull of the polytope dened by generation limits and the minimum up/downtime constraints that will be presented shortly. With respect to constraints (2.1d), the use of p 0 g;t variables necessitate an update in the demand constraints, as below: X g2G p 0 g;t + C g x g;t d t 8t2T: (2.5) Finally, the following inequalities ensure that the reserve requirements can be fullled: X g2G r g;t t 8t2T: (2.6) CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 12 Formulating operational and technological restrictions In order to compute realistic production schedules, the MIP formulations must take into account the physical limitations of the components of the power systems. Among these limitations, some of the most essential ones are the minimum up/downtime and ramping restrictions. The minimum up/downtime restrictions ensure that the on/o status of the generators do not change rapidly. Frequent state transitions have several adverse consequences including (i) increased operator stress, (ii) diminished generator life, and (iii) increased emission of pollutants during transient periods (see Takriti and Birge (2000)). Such restrictions are quite practical and are included in many commercial tools. To formulate these restrictions, several constraints were proposed in Arroyo and Conejo (2000); Takriti and Birge (2000); Lee et al. (2004). Following Rajan and Takriti (2005), the minimum uptime and downtime constraints are best formulated using the following constraints: t X i=tUTg +1 s g;i x g;t g2G; t2T; (2.7a) t X i=tDTg +1 s g;i 1x g;tDTg g2G; t2T: (2.7b) Observe that when the generator g is operational at time t, the right-hand side of constraint (2.7a) is set to 1. In this case, the generator may have been turned on at most once in the last UT g periods (due to minimum uptime restrictions). On the other hand, if it is idle at timet, it could not have been turned on in the lastUT g time periods, as otherwise it should be operational at timet. Constraint (2.7b) is just a rewritten version of a similar constraint for the downtime requirements. It has been widely observed that these turn on/o inequalities signicantly outperform its contenders Ostrowski et al. (2012). Indeed, Rajan and Takriti (2005) showed that constraints (2.7) dene facets of the polytope dened by the minimum up/downtime constraints. The ramping restrictions limit the maximum change in production between consecutive periods and ensure that the generation requirements can be matched by the electricity production without exceeding the generator limitations over extended periods of time. A basic representation of these ramping restrictions, in the space of p 0 g;t and r g;t variables, appears below: p 0 g;t +r g;t p 0 g;t1 R g 8g2G; t2T; (2.8a) p 0 g;t1 p 0 g;t R g 8g2G; t2T: (2.8b) The above inequalities, by themselves, do not take into account the start up and shut down rates of the generators. Indeed, certain generation limit constraints (such as (2.4)) must be used in conjunction with (2.8) to make sure these restrictions are also satised. The ramping inequalities may also be written in the space of p g;t and p g;t variables. Below, we give the formulation of ramping restrictions that appeared in Ostrowski et al. (2012). p g;t p g;t1 S g s g;t + R g x g;t1 8g2G; t2T; (2.9a) p g;t1 p g;t S g z g;t + R g x g;t 8g2G; t2T: (2.9b) Constraint (2.9a) ensures that generator g cannot ramp up more than S g if it has just been turned on, or R g if it remains on at time t. Similarly, constraint (2.9b) limits the decrease in the power output by R g at any time that the generator remains operational. If the plant is turned o att (x g;t = 0), then the output of the generator cannot be larger than S g to obey the shut down limits. Observe that (2.9a) cannot be tight if the generator is turned o at timet ( p g;t = 0) becausep g;t1 C g < 0, but the right-hand-side is R g > 0. A similar argument also applies to (2.9b) when the generator is turned on at time t. In fact, these (along with many others) were the motivation behind studying tighter ramping constraints and valid inequalities in Ostrowski et al. (2012); Damc-Kurt et al. (2015); Jiang et al. (2016); Pan and Guan (2016). CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 13 In particular, Damc-Kurt et al. (2015) provided the following strengthened ramping constraints: p g;t p g;t1 ( S g R g C g )s g;t + ( R g + C g )x g;t C g x g;t1 8g2G; t2T; (2.10a) p g;t1 p g;t ( S g R g C g )z g;t + ( R g + C g )x g;t1 C g x g;t 8g2G; t2T: (2.10b) These inequalities were proved to be facet-dening for the two-period ramp-up and ramp-down polytopes (i.e., the polytopes of the UC problem that are limited to two consecutive periods and consider only the ramp-up and ramp-down constraints, respectively). Linearization of the objective function The nonlinear operating cost functions in the objective of a UC problem are typically approximated with piecewise-linear convex functions. We begin with the linearization of the start up costs. For a generator g, we use the setI g = (t 1 g :::t jIgj g )T to denote the set of idle periods after which the cost incurred to turn on a generator changes. The start up costs are given by fc g where 2I g . These costs are assumed to obey fc jIgj g :::fc 1 g ,8g2G, indicating that the start up costs tend to increase as the idle time of the generators grow. Accordingly, a generator that has been idle for 1 to t 1 g 1 periods will incur a start up cost of fc 1 g , t 2 g to t 3 g 1 periods will incur fc 2 g , and so forth. It is easy to see that the start up costs are determined by a step-function. This function can be incorporated into an MIP formulation using the approaches in Nowak and R omisch (2000) and Morales-Espa~ na et al. (2013), where the former linearizes it and the latter partitions its domain. These approaches respectively require the following variables: f g;t : Incurred start up cost for generator g at time t, g;t; : Start up cost selection variable (1 if g must incur the start up cost fc g at time t, 0 otherwise). Using thef g;t variables, Nowak and R omisch (2000) determine the cost of turning on a generator with the following start up cost constraints: f g;t fc g x g;t t g X i=1 x g;ti 8g2G;t2T;2I g : (2.11) To minimize the total start up costs, the objective must then contain the sum of all f g;t variables. Alternatively, the approach in Morales-Espa~ na et al. (2013) formulates the start up costs as follows: g;t; t +1 g 1 X i=t g z ti 8g2G; t2 t +1 g ;:::;jTj ; 2I g n 1:::DT g 1; t jIgj g ; (2.12a) jIgj X =1 g;t; =s g;t 8g2G; t2T: (2.12b) Constraints (2.12) ensure that a single and correct start up type is selected based on how long the unit has been idle. Notice that (2.12a) need not be dened for < DT g (since an idle generator must obey the minimum downtime restrictions) and for =t jIgj g (as it will be redundant due to (2.12b)). The integrality restrictions on g;t; need not be enforced provided that the start up costs are monotonically increasing with the number of periods the generator remains idle. In general, the production costs are formulated with piecewise-linear convex functions, as the marginal cost of production increases with increasing levels of production. We dene a set of production amounts p g (with 2f1::: max g g) at which the incurred unit cost of production for generatorg changes. These production levels can be interpreted as the breakpoints where the slope of the piecewise-linear function is altered. The unit generation cost within the interval p g ;p +1 g is denoted withvc g . Similar to the start up costs, we assume thatvc max g g ::: vc 1 g ; 8g2G, which suggests that the marginal costs increase as the outputs approach the generator capacities. Finally, the aggregate cost of generating p g units of output is represented through the function ^ V (p g ). Following these denitions, the production amounts can be accounted by the variables v g;t with the following production-cost CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 14 constraints: v g;t vc g (p g;t p 1 g ) + ^ V g (p 1 g ) 8g2G; t2T; and8 =f1::: max g g: (2.13) Above, we assume that p 0 g = 0 (and hence, ^ V g (p 0 g ) = 0). Along with its other benets, the use of p 0 g;t variables also eases the formulation of piecewise-linear production costs. Recall that the total production of a generator is now accounted by two terms,p 0 g;t and C g x g;t . As long as the generator is operational, the cost of producing the initial C g units of energy is xed and can be computed a-priori (herein denoted as ^ V g ( C g )). This cost can be associated with the s g;t and ~ x g;t variables and directly accounted in the objective function. Therefore, the variables v g;t now provides the cost of producing an amount of energy in addition to C g . This cost is computed by the following production-cost constraints: v g;t vc g (p 0 g;t + C g p 1 g ) + ^ V g (p 1 g ) ^ V g ( C g ) 8g2G; t2T; and8 =fl g ::: max g g: (2.14) Above, l g ( 1) is the value of the index for which p lg1 g = C g . For values of 2f1:::l g 1g, constraints (2.14) are omitted because their right-hand side will be no larger than zero and will be trivially satised due to the nonnegativity of the variables. 2.1.2. Benchmark Formulations Using the inequalities presented in the previous section, we provide two complete UC formulations, which are primarily based on Ostrowski et al. (2012) and Morales-Espa~ na et al. (2013), respectively. The formulations are named after the initials of their corresponding authors. OAV MLR Description min P t2T P g2G f g;t +v g;t P t2T P g2G C g x g;t +v g;t objective functions, + P 2Ig fc g g;t; subject to: (2.1a) (2.1a) state constraints; (2.1b); (2.2) (2.4) production limits; (2.1d); (2.3) (2.5); (2.6) demand and reserve requirements; (2.7) (2.7) min up/downtime requirements; (2.9) (2.8) ramping restrictions; (2.11); (2.13) (2.12); (2.14) start up and production costs; (x; s; z)2f0; 1g 3jGjjTj (x; s; z)2f0; 1g 3jGjjTj integrality restrictions, (p; p; f; v) 0 (p 0 ; r; ; v) 0 nonnegativity restrictions. 2.2. The State-Transition Formulation In this section, we present the state-transition formulation (STF) for the UC problem. We develop this formu- lation through the use of the following state-transition variables: ~ x g;t : 1 if g remains operational at time t, 0 otherwise. ~ s g;t =~ z g;t : 1 if g is turned on / o at time t, 0 otherwise, To see how these variables capture the transition of generator states, we provide an illustration in Figure 2.1. In this gure, the on/o status of a generator is denoted with nodes and the feasible state transitions are represented with arcs. Notice that the variables for the decision remaining o has not been dened and the corresponding arc in Figure 2.1 is marked with a dashed line. This is because such a transition is completely determined by the values of the other state-transition variables. Indeed, exactly one state transition occurs at each time period, and the corresponding decision variable is set to 1. If none of these variables are set to 1, then the remaining o transition CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 15 corresponding to the dashed line is said to occur. Following this discussion, it is easy to see that the inequality ~ s g;t + ~ x g;t + ~ z g;t 1 8g2G; t2T: (2.15) is true for any feasible generator schedule formulated with state-transition variables. on t-1 off t-1 on t off t ˜ x t ˜ s t ˜ z t Figure 2.1: Illustration of the state-transition variables. A deeper look at Figure 2.1 reveals the network sub-structure embedded in the new formulation. In particular, we observe that a sequence of state transitions of a single generator can be de facto represented by a path on the graph comprising states as nodes and state transitions as arcs. Writing out the ow conservation constraints at the nodes of this graph, we immediately obtain the following state-transition constraints for the UC problem: ~ s g;t1 + ~ x g;t1 = ~ z g;t + ~ x g;t 8g2G; t2T: (2.16) Constraints (2.16) are the ow conservation constraints at the nodes where the state of the generator is on. At the o nodes, the ow conservation constraints are the same, therefore omitted. This can be easily veried by assigning the expression 1 ~ s g;t ~ x g;t z g;t to the dashed arc in Figure 2.1. From our perspective, these constraints indicate that if a generator is operational at time t 1 then it must either be turned o or remain on at time t. The equation (x g;t ) = (~ s g;t + ~ x g;t ) allows us to translate any constraint derived for the benchmark formulations into the constraints of our formulation. Although the constraints dier in essence, it is easy to verify that (2.16) can also be derived from (2.1a) using the transformation described above. In a UC problem, the above formulation of the state transitions can potentially lead to signicant computational gains, especially when large numbers of state transitions are anticipated. Accordingly, for problems involving longer planning horizons and facing variable demand patterns, one should expect a better performance from this new formulation, compared to contenders that use (x; s; z) variables. While other components of a UC problem can be formulated in alternative ways, the constraints presented in the rest of this section perform well in tandem, and some naturally dene the facets of their corresponding polytopes. Formulating production In our formulation, we opt for the p 0 g;t and p g;t variables for formulating production amounts and reserve requirements. With the state-transition variables, the on/o status of a generator is given by the expression ~ s g;t + ~ x g;t . Accordingly, the production limits for each g2G; t2T can be given as follows: p g;t p 0 g;t + C g (~ s g;t + ~ x g;t ); (2.17a) p g;t C g (~ s g;t + ~ x g;t ) + ( S g C g )~ z g;t+1 : (2.17b) Constraints (2.17a) ensure that the maximum possible generation amount at time t ( p g;t ) is greater than the actual production amount at the same time period, and constraints (2.17b) bound p g;t by the total capacity or the shut down limit of the generator. The minimum production constraints are redundant due to the use of thep 0 g;t variables. As is the case in (2.2), the variable p g;t must also respect the start up limits, but this restriction is omitted in (2.17) as it will be implied by the ramping constraints. CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 16 The demand constraints of the STF are similar to that of MLR given in (2.5): X g2G p 0 g;t + C g (~ s g;t + ~ x g;t ) d t 8t2T: (2.18) Finally, the reserve requirements are formulated without eort by adopting the constraints in (2.3). Formulating operational and technological restrictions We now formulate the minimum up/downtime and the ramping restrictions. We begin with the former, and translate the minimum up/downtime inequalities (2.7) into our formulation as presented below: t1 X i=tUTg +1 ~ s g;i ~ x g;t 8g2G; t2T; (2.19a) t X i=tDTg ~ s g;i 1 ~ x g;tDTg 8g2G; t2T: (2.19b) Constraint (2.19a) ensures that if the generator remains on, it could have been turned on at most once in the previous UT g 1 time periods. If it does not remain on, then it could not have been turned on in these time periods due to minimum uptime restrictions. Similarly, constraint (2.19b) ensures that if the generator remains on, it cannot be restarted in the current time period or in the next DT g time periods due to minimum downtime restrictions. On the other hand, if it does not remain on, it can be turned on at most once in these time periods, due to minimum uptime and downtime restrictions. To formulate the ramping restrictions, we take advantage of both the state-transition and the production vari- ables, and provide the following ramp up constraints: p g;t p 0 g;t1 S g ~ s g;t + ( R g + C g )~ x g;t 8g2G; t2T: (2.20a) The coecient of the remain-on variable is increased by C g because thep 0 g;t1 variable only accounts for the power generation beyond C g . When a generator is turned on in periodt, it should have been idle in periodt1. Therefore p 0 g;t1 = 0 and no increment for the ~ s g;t coecient is necessary. The ramp down constraints are formed in a similar manner and given below. p 0 g;t1 p 0 g;t ( S g C g )~ z g;t + R g ~ x g;t 8g2G; t2T: (2.20b) Whenever the generator remains on, the above inequality limits the reduction in the power output such that the ramp down restrictions are obeyed. When the generator is turned o at time t, the inequality reduces to a tight upper bound on the p 0 g;t1 variable. When the reserve requirements are ignored, it can be shown that (2.20) is equivalent to (2.10). Linearization of the objective function For the linearization of the start up costs, we opted for the use of the f g;t variables. Observe that, when turned on, the start up cost of a generator will at least be its warm-start cost, i.e., the start up cost incurred when the generator has not cooled down since the previous operational state. Due to the minimum downtime restriction, this cost is at least fc DTg g . Furthermore, it can be associated with the start up variable and accounted directly in the objective function using the additional term fc DTg g ~ s g;t . Therefore, the variablef g;t will now represent the extra cost of turning on a generator which has been idle for some time that is longer than its minimum downtime (DT g ). In view of this observation, variables f g;t must satisfy the following start up cost constraints: f g;t (fc g fc DTg g ) ~ s g;t X i=DTg ~ s g;ti ~ x g;t 8g2G; t2T; and 2I g nf1:::DT g g (2.21) CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 17 We note that, in general, the number of constraints to formulate the piecewise-linear cost functions can make mixed-integer programs extremely large. As in MLR, the representation above mitigates this potential issue by eliminating up tojTj P g2G DT g constraints. Finally, the production costs can be formulated by directly incorporating (2.14) into our formulation. Similar to MLR, the cost of minimum required production is determined by C g (~ s g;t + ~ x g;t ), and appended to the objective function. The state-transition formulation. A summary of the new formulation is provided below: min P t2T P g2G f g;t +v g;t + (fc DTg g + C g )~ s g;t + C g ~ x g;t subject to: (2.16) state transition constraints; (2.17) production limits; (2.3); (2.18) demand and reserve requirements; (2.19) min up/downtime constraints; (2.20) ramping restrictions; (2.21); (2.14) start up and production cost functions; (~ x; ~ s; ~ z)2f0; 1g 3jGjjTj integrality restrictions; (p; p; f; v) 0 nonnegativity restrictions. In STF, we do not include constraints (2.15) as they are implied by other components of our formulation. A formal argument is provided in the next section. At this point, we nd it convenient to summarize all the notation used in §2.1 and §2.2 in Table 2.1. Sets G: Set of generators. T : Set of discretized time periods. Parameters dt: Demand in period t2T . t: Reserve requirements in period t. Cg : Production capacity of g2G. Cg : Minimum required production amount wheng is in operational state. Rg : Ramp up limit of g. Rg : Ramp down limit of g. Sg : Ramp up limit of g when g is turned on. Sg : Ramp down limit of g before g is turned o. UTg : Minimum number of periods that g must remain operational after it has been turned on. DTg : Minimum number of periods that g must remain idle after it has been turned o. Ig : Set of idle-time durations after which the start up cost of g increases. fc g : Start up cost of g when g remains idle for 2Ig periods. p g : Production amount after which the unit cost of pro- duction increases (2f1::: max g g). vc g : Unit cost of production when the production amount lies within [p g ;p +1 g ). ^ V (p g ): Total cost of producing p g units of output by g. Decision Variables xg;t: State variable (1 ifg is operational at timet, 0 oth- erwise). sg;t: Start up variable (1 if g is turned on at time t, 0 otherwise). zg;t: Shut down variable (1 if g is turned o at time t, 0 otherwise). ~ xg;t: Remain-on variable of STF (1 if g remains opera- tional at time t, 0 otherwise). ~ sg;t: Start up variable of STF (1 ifg is turned on at time t, 0 otherwise). ~ zg;t: Shut down variable of STF (1 if g is turned o at time t, 0 otherwise). pg;t: Amount of production by g at time t. p 0 g;t : Variable amount of production by g at time t (ex- cludes Cg ). pg;t: Production amount thatg can supply at timet (in- cludes the amount committed to fullling reserve re- quirements). rg;t: Production amount that g can supply for fullling reserve requirements. fg;t: Start up cost that g incurs at time t. g;t; : Start up cost selection variable (1 if g must incur start up cost fc g , 0 otherwise). vg;t: Production cost that g incurs at time t. Useful Relations pg;t =pg;t +rg;t pg;t = Cgxg +p 0 g;t xg;t = ~ xg;t + ~ sg;t Table 2.1: Nomenclature for the mathematical formulations STF, OAV and MLR. CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 18 2.2.1. Insights on the Polyhedral Structure This section summarizes the polyhedral features of our formulation and draw connections with other formulations from the literature. Our discussion will focus on modeling components that involve a single generator, therefore the generator indices are neglected for notational ease. We begin with showing the redundancy of (2.15) in STF. Proposition 2.1. Constraints (2.15) are implied by the minimum downtime constraints (2.19b) and the state- transition constraints (2.16). Proof. Consider constraint (2.19b) for time t +DT 1, 1 ~ x t1 + t+DT1 X i=t1 ~ s i = ~ x t1 + ~ s t1 + t+DT1 X i=t ~ s i [by (2.16)] = ~ z t + ~ x t + ~ s t + t+DT1 X i=t+1 ~ s i ~ z t + ~ x t + ~ s t : Remark 2.1. If DT = 1, constraints (2.15) become equivalent to the minimum downtime constraints. The following denitions and facts will be useful in the remainder of this section. We use dim(P) to denote the dimension of a polyhedronP. A facetF ofP has dim(F) = dim(P) 1. Accordingly, to prove that a constraint | x 0 denes a facet ofP, it is sucient to identify dim(P) anely-independent points (x i ) which satisfy | x i = 0 ;8i = 1::: dim(P). A facet-dening inequality will dominate all other constraints, therefore sought in most polyhedral analyses. For proofs of this type appealing to ane independence for facet-dening inequalities, the reader may refer to Nemhauser and Wolsey (1999). Next, we show that the use of the state-transition decision variables does not aect the polyhedral properties of the original minimum up/downtime constraints (2.19) of Rajan and Takriti (2005), proposed for formulations with (x; s; z) variables. Consider the minimum up/downtime polytope of the STF: M = ~ s; ~ x;~ z2f0; 1g jTj : P t1 i=tUT+1 ~ s i ~ x t 8t2fUT:::jTjg, P t i=tDT ~ s i 1 ~ x tDT 8t2fDT + 1:::jTjg, ~ s t1 + ~ x t1 = ~ z t + ~ x t 8t2f2:::jTjg: In the denition above, the history of the generators is neglected by removing the corresponding turn on/o constraints dened for periodsf1:::UT 1g, andf1:::DTg, respectively. Proposition 2.2. The minimum up/downtime constraints (2.19) are the facets of the polytope conv(M). Proof. Observing the linear independence of equations in (2.16), dim(conv(M)) can at most be 2jTj + 1. Indeed this is true because the following 2jTj + 2 integer solutions ofM are anely independent. To see this, note that by ordering the solution in a matrix and placing rst the ~ x variables and then the ~ s and the ~ z variables, we obtain a lower-triangular submatrix. Consider a set ofjTj solutions (~ s; ~ x;~ z) i , where i2f1:::jTjg. This solution has the following structure. The components of the vector s (resp. z) are set to zero with the only exception of thei th component (resp. (i+UT +1) th orjTj, whichever is smaller). This element is set to 1. The components of ~ x are set as ~ x k = ( 1; i + 1k minfi +UT;jTjg 0; otherwise: CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 19 jTj additional solutions (~ s; ~ x;~ z) j ; j2f1:::jTjg, which belong to a set J, are obtained by setting the ~ s and ~ z vectors to null, while the vector ~ x has the following components: ~ x k = ( 1; kj 0; otherwise: The last two solutions are (~ s; ~ x;~ z) = (0; 0; (1; 0;::: 0)) and the null solution. Given the dimension of conv(M), any constraint (2.19a) is a facet if it has dimension 2jTj (i.e., if there are 2jTj + 1 anely independent solutions which satisfy the inequality as equality). For some t2fUT:::jTjg, consider the set U t = n ~ s; ~ x;~ z2f0; 1g jTj : P t1 i=tUT+1 ~ s i = ~ x t o : All the solutions listed above belong to the setU t , with the exception ofjTj (t + 1) solutions of set J. In these solutions, the generator is turned o at time t + 1 or later (i.e., all solutions (~ s; ~ x;~ z) j where j > t). For each of these solutions, turn the generator on at time tUT + 1 by setting the (tUT + 1) th component of the vector s to 1, and setting the rst tUT + 1 components of the vector ~ x to 0. 2jTj + 1 anely independent solutions have thus been generated, proving thatU t is a facet of conv(M). The proof for minimum downtime constraints (2.19b) is analogous and can be obtained with the same arguments made for constraints (2.19a). Next, we show the validity of the ramping inequalities (2.20), and demonstrate their relation to their counterparts in Ostrowski et al. (2012). Proposition 2.3. For each t2T , the ramp-up inequality (2.20a) is valid and dominates constraint p t p t1 S~ s t + R~ x t ; (P1a) which is equivalent to (2.9a) when the state-transition variables are used. Proof. It is straightforward to verify the validity of constraints (2.20a). The inequality limits the increase in production by S if the generator has just been turned on (in this case p 0 t1 = 0), or by ( R + C) if it remains on. As also pointed out in §2.1.1, the ramp-up inequality (P1a) is inactive if the generator has just been turned o; but it can be strengthened by lifting it to the space of ~ z t variables. p t h p 0 t1 + C(~ s t1 + ~ x t1 ) i S~ s t + R~ x t C~ z t : (P1b) The term within brackets is p t1 . To verify that (P1b) is valid, it is sucient to consider the case where ~ z t = 1 (the case ~ z t = 0 is trivial). As the generator is not operational at time period t, p t is 0 and the inequality reduces top 0 t1 + C(~ s t1 + ~ x t1 ) =p t1 C, which is the minimum production requirement. Therefore (P1b) is valid and dominates (P1a), which is easy to see as the left-hand side of (P1a) assumes no smaller value than the left-hand side of (P1b). Constraints (2.20a) can then be obtained from (P1a) using (2.16): p t p 0 t1 S~ s t + R~ x t C~ z t + C(~ s t1 + ~ x t1 ) = S~ s t + R~ x t C~ z t + C(~ z t + ~ x t ) = S~ s t + ( R + C)~ x t : The analysis for the ramp down constraints is provided next. Proposition 2.4. For each t2T , the ramp-down inequality (2.20b) is valid and dominates the seed inequality p t1 p t S~ z t + R~ x t , which is equivalent to (2.9b) when the state-transition variables are are used. CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 20 Proof. Consider the seed inequality: h p 0 t1 + C(~ s t1 + ~ x t1 ) i h p 0 t + C(~ s t + ~ x t ) i S~ z t + R~ x t ; which can be strengthened by a lifting to the space of ~ s t , h p 0 t1 + C(~ s t1 + ~ x t1 ) i h p 0 t + C(~ s t + ~ x t ) i S~ z t + R~ x t C~ s t : When ~ s t = 1, the power output at timet1 must be 0, and the rst term in brackets disappears. Since ~ x t = ~ z t = 0, the inequality simply reduces to p 0 t 0, which conrms its validity. By simple algebra, we obtain p 0 t1 p 0 t S~ z t + ( R + C)~ x t C(~ s t1 + ~ x t1 ) = S~ z t + ( R + C)~ x t C(~ z t + ~ x t ) (P2) = ( S C)~ z t + R~ x t where (P2) is achieved using (2.16). Using the relations (2.1a) and (2.16), it is trivial to verify the following remark. Remark 2.2. When the reserve requirements are neglected, the ramping inequalities (2.20a) and (2.20b) of the new formulation are equivalent to the facet-dening two-period ramping inequalities of the base formulation, proposed by Damc-Kurt et al. (2015), i.e., constraints (2.10a) and (2.10b). In contrast, the ramping constraints in MLR are such that the state variables do not explicitly guide the ramping rates of generators. For our nal analysis, we introduce the two-period ramping polytope in our formulation, under the assumption thatUT > 1. For clarity of exposition, we rst give the generation and ramping constraints for a given t as follows: p t p 0 t + C(~ s t + ~ x t ) (2.23a) p t S~ s t + C~ x t + ( S C)~ z t+1 (2.23b) p 0 t1 p 0 t ( S C)~ z t + R~ x t (2.23c) p t p 0 t1 S~ s t + ( R + C)~ x t + ( S R C)~ z t+1 : (2.23d) It is easy to see that (2.23b) is a tightened version of (2.17b), whereas (2.23d) is the lifted version of the ramp up inequality (2.20a). Our computations do not consider (2.23b) and (2.23d) in order to simplify the exposition, and avoid detracting from the paper's main messages. It should also be noted that (2.23d) is equivalent to the valid inequality (23) of Ostrowski et al. (2012), when written in terms of the state-transition and p 0 g;t1 variables. We give the two-period ramping polytope below. R 2 t = n ( p t ; p 0 t1 ; p 0 t ; ~ s t ; ~ x t ; ~ z t ~ z t+1 )2R 3 + f0; 1g 4 : (2.23a); (2.23b); (2.23c); (2.23d); ~ s t + ~ x t + ~ z t 1; ~ z t+1 ~ x t o With respect to the two-period ramping polytope dened in Damc-Kurt et al. (2015), the above denition has some dierences. For instance, we consider both the ramp up and ramp down constraints together, which were treated separately in Damc-Kurt et al. (2015). Furthermore, our polytope is dened in the space of p t , p 0 t , p 0 t1 , ~ s t , ~ x t , ~ z t and ~ z t+1 variables, and considers the reserve requirements. Proposition 2.5. Assuming C > R + R 2 R and the generation limit constraints, (2.23a) and (2.23b), and the ramping down inequality (2.23c) and (2.23d) are facets of the two-period ramping polytope conv(R 2 t ). Proof. The dimension of conv(R 2 t ) can at most be seven. In Table 2.2, we present feasible solutions for this polytope. The rst eight solutions are anely-independent, thus proving the full dimensionality of conv(R 2 t ). Using the solutions displayed in Table 2.2, it is also easy to verify that constraints (2.23a), (2.23b), (2.23c) and (2.23d) are facet-dening for conv(R 2 t ). In the last column of Table 2.2, we report the constraints which are CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 21 Table 2.2: Set of feasible solutions for the two-period ramping polytope. # pt p 0 t1 p 0 t ~ st ~ xt ~ zt ~ z t+1 Inactive Cons. 1 0 0 0 0 0 0 0 2 0 0 0 0 0 1 0 (2.23c) 3 0 S C 0 0 0 1 0 (2.23d) 4 S 0 0 1 0 0 0 (2.23a) 5 S 0 S C 1 0 0 0 (2.23c) 6 C C C R 0 1 0 0 (2.23a); (2.23d) 7 R + C 2 R R 0 1 0 0 (2.23b); (2.23d) 8 S S + R C S C 0 1 0 1 (2.23d) 9 R + C R 0 0 1 0 0 not active in the corresponding solution. To prove that each of the constraints are facet-dening, we modify the coecients of the solutions in Table 2.2 as follows. • Constraint (2.23a): change the coecient of p t from S to C in Solution 4. • Constraint (2.23b): seven solutions out of the eight already satisfy the constraint as equality. • Constraint (2.23c): substitute Solution 5 with Solution 9. • Constraint (2.23d): change the coecients of p 0 t1 from C to C R C, from 2 R to 0, and from S + R C to S R C in Solutions 6, 7 and 8 respectively. For each of the constraints, we have provided six anely-independent feasible solutions which satisfy all constraints as equality, thus proving the proposition. Remark 2.3. Observe that for every t2T , the inequality (2.15) is also facet-dening for the polytope conv(R 2 t ). 2.3. Computational Experiments Our analysis will focus on two classes of instances. The rst stems from the synthetic instances developed by Ostrowski et al. (2012), which are based on Carri on et al. (2006). These are single-day instances, involving 24 time periods, and are characterized by increasing numbers of generators. They are known to be challenging for branch- and-bound algorithms, and therefore, have been commonly used in the literature for computational experiments. Using these instances, we generated two additional data sets. In the rst, the numbers of generators are increased tenfold, and in the latter, the time horizon is extended to seven days. Instances with a longer time horizon are more challenging not simply because of the increased dimensions but also due to the daily trends and uctuations in the demand. The second class of instances are based on the realistic instances obtained through the FERC. Further details can be found in Krall et al. (2012). For both classes of instances, details are given in Appendix A and the references therein. For dynamic planning models it is well known that the planning horizon can have a signicant impact on the decisions. The resulting decisions are often myopic, and this phenomenon is commonly referred to as the planning horizon eect. Such myopic choices can be remedied by choosing to solve longer-horizon models, and only implementing the day-ahead plan. This approach is sometimes also referred to as a receding (rolling) horizon approach (RHA). In case of the UC model, since a week can be considered as a regeneration point for demands (see Dupa cov a et al. (2003); Sen et al. (2006)), such an RHA is likely to avoid myopic choices. In keeping with this outlook, the choice of a 168 periods re ects the UC problem for the RHA. In Table 2.3, we present statistics regarding the considered formulations. We observe that the number of binary variables is the same across all formulations. In terms of the numbers of variables, constraints, and nonzero coecients, STF attains the minimum amounts, leading to compact descriptions of UC problems. We note that the week-long realistic instance contains a signicantly larger set of generators compared to other week-long instances that we have experimented with. CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 22 Synthetic Instances (averaged) Realistic Instances jTj Binary / Total Vars. Constraints Nonzeros Binary / Total Vars. Constraints Nonzeros OAV 24 77,076 / 179,845 304,789 1,178,779 67,608 / 157,753 287,108 1,050,268 168 53,953 / 125,892 214,676 868,319 473,256 / 1,104,265 2,016,548 9,327,706 STF 24 77,076 / 179,845 231,313 1,012,799 67,608 / 157,753 225,188 925,739 168 53,953 / 125,892 163,243 748,272 473,256 / 1,104,265 1,583,108 8,416,361 MLR 24 77,076 / 205,537 257,125 1,048,143 67,608 / 180,289 241,532 945,153 168 53,953 / 143,876 181,312 769,461 473,256 / 1,262,017 1,697,516 8,519,871 Table 2.3: Formulation statistics (synthetic instances withjTj = 24 contain ten times more generators) All runs were performed on a single thread of a Dell Desktop PC with Intel ® Core ™ i7-3770S CPU @ 3.10 GHz, 7.68 GB of RAM, and running Ubuntu Linux 12.04.3 LTS. The formulations were solved with CPLEX 12.5.1. The default parameters of CPLEX were preserved, but a time limit of 2 hours is imposed. For instances which could not be solved within this limit, we report the relative optimality gap based on the best available solution and lower bound. It is important to note that the benchmark formulations are true to the mathematical representations provided by the original authors. In making comparisons, certain choices -such as software parameters- are kept the same across all formulations, and no tuning was performed to individually improve their performances. In Tables 2.4 and 2.5, we report the solution times for the synthetic instances, along with the numbers of cuts and nodes generated by the solver. When accompanied with fast solution times, small numbers of cuts could indicate that the instances are amenable to producing integer solutions fast and with little need for tightening. Likewise, small numbers of branch-and-bound nodes (or the lack thereof) imply that the non-convexity in the instances can be easily tamed. We rst make comparisons on the single-day instances. We observe that STF and MLR perform considerably faster than OAV, especially as the number of generators grow. In terms of solution times, there does not appear to be signicant dierences among STF and MLR. However, it is promising to see that the total solution time for all instances is the smallest when STF is used, and the individual times never exceed a minute. Moreover, with STF, all of the instances were solved at the root node, with the minimum number for cutting planes. In contrast, we observe that multiple nodes were explored in a few instances using the benchmark formulations. OAV STF MLR jGj Cuts Nodes Time Cuts Nodes Time Cuts Nodes Time 280 - - 2.5 - - 1.9 - - 1.1 350 506 - 5.2 - - 2.7 - - 1.6 440 3,729 - 11.0 - - 3.5 - - 2.2 450 3,867 - 34.4 218 - 5.0 1,019 100 12.4 490 6,334 - 29.8 - - 4.7 63 - 5.9 500 9,520 - 83.5 863 - 8.5 2,142 80 24.7 510 1,850 - 18.0 - - 4.3 2 - 4.6 510 6,080 - 23.5 - - 4.1 - - 2.8 520 6,043 - 26.8 - - 4.6 1 - 3.8 540 7,299 - 36.0 - - 4.6 12 - 4.4 1320 2,676 - 35.9 - - 16.1 - - 9.8 1560 12,059 - 187.2 - - 23.4 1 10 33.4 1560 12,737 - 111.6 - - 21.6 - - 19.4 1650 17,249 - 133.1 - - 28.6 - - 23.9 1670 19,588 - 230.5 - - 29.6 - - 27.0 1720 7,642 - 81.3 - - 26.4 - - 18.9 1820 22,170 524 1,069.9 - - 33.8 753 10 61.9 1820 24,087 250 944.3 36 - 49.1 306 485 263.7 1830 15,978 - 281.3 - - 36.7 2 - 41.8 1870 24,373 - 508.8 1,386 - 58.7 3,636 - 65.9 Total 203,787 774 3,854.3 2,503 - 367.8 7,937 685 629.2 Table 2.4: Performance measures of branch-and-bound algorithm for synthetic instances (jTj = 24). When we consider week-long instances, we observe an important feature of our formulation. Recall from §2.2 that the STF involves an implicit network sub-structure, and it has better potential in capturing the state transitions CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 23 OAV STF MLR jGj Cuts Nodes Time Cuts Nodes Time Cuts Nodes Time 28 428 - 3.3 12 - 3.4 252 170 5.8 35 914 - 6.4 - - 2.1 865 - 4.4 44 1,188 - 7.9 - - 2.8 684 - 4.9 45 3,016 417 34.3 158 - 4.7 984 413 12.5 49 2,981 210 25.1 92 75 8.7 517 506 21.3 50 3,191 529 61.2 216 - 6.0 837 806 23.1 51 1,263 - 11.3 18 - 5.1 105 - 4.9 51 1,550 - 11.0 - - 3.1 273 185 8.8 52 2,099 - 14.5 100 - 4.7 469 - 6.4 54 2,470 - 14.7 88 63 8.2 373 739 21.9 132 1,206 - 32.9 - - 11.9 2,159 - 18.8 156 2,781 - 60.7 51 - 23.5 285 200 39.6 156 3,413 - 48.3 - - 15.7 1,994 342 52.0 165 6,564 - 102.8 - - 21.3 3,306 531 89.6 167 8,292 - 111.1 243 - 22.1 638 1,880 432.2 172 2,705 - 75.4 310 - 33.5 2,522 - 28.2 182 7,556 - 103.6 189 - 24.6 2,469 150 51.8 182 11,135 - 137.6 629 - 24.9 5,404 624 76.9 183 9,038 - 129.7 285 - 24.4 1,887 577 110.3 187 10,979 140 205.9 389 - 35.4 3,209 120 60.4 Total 82,769 1,296 1,231.8 2,591 138 286.3 26,763 7,093 1,098.8 Table 2.5: Performance measures of branch-and-bound algorithm for synthetic instances (jTj = 168). of the generators between consecutive time periods. Aligned with our expectations in §2.2, we observe that the complexity introduced by a longer time horizon is best tamed with STF. The solution times are signicantly better than the benchmark formulations for almost every instance. Furthermore, there are only two instances for which multiple branch-and-bound nodes have been explored. In contrast, the benchmark formulations needed to explore hundreds of nodes to optimize a signicant portion of the instances. These results reveal the potential of our formulation to scale up to more realistic instances solved by the power industry. We turn our attention to our realistic instances. We consider both a day-long and a week-long version of the instance made available by Krall et al. (2012). The analysis of these instances will provide better intuition on how the new and the benchmark formulations perform in the problems solved in the energy industry. Comparing Table 2.6 with Tables 2.4 and 2.5, we observe that the dierences in the computational performance of the new model and the benchmark formulations are more pronounced. For instance, we observe that OAV hits the time limit on the week-long problem, whereas MLR spends more than 20 minutes to obtain an optimal solution. In comparison, the new model can provide an optimal solution within 12 minutes, achieving a 45% reduction. We conclude our analysis by discussing the percentage of fractional variables in Table 2.6. A small number of fractional variables in the continuous relaxation of the problem may serve as an indication of the tightness of a formulation. Observe that the new formulation provides far larger numbers of integer variables in the LP-relaxations of the problems, conrming the eectiveness of the formulation. The objective value of the LP-relaxation of STF can be 1.4% and 0.2% better than that of OAV and MLR, respectively. More importantly, in the root relaxation (i.e., the starting LP-relaxation created within CPLEX), the percentage of fractional variables shrink even further. Indeed, in the week-long instance, 99.82% of the integer variables are observed to be integral, leading to an almost-integer (lower-bounding) solution. In order to further assess the behavior of the formulations, we have experimented with modied versions of the original synthetic and realistic instances. In particular, the hourly demand and reserve data are multiplied with (1+ ), where is chosen to mimic realistic changes in demand (see New York ISO (2016)). Likewise, generator capacities in synthetic instances are perturbed with the same coecient (1 +). For realistic instances, however, generator capacities are not altered in order to preserve their authenticities. Table 2.7 gives a summary of our analysis (where solution times are summed for all 20 synthetic instances). Examining the results for synthetic instances, we observe that STF is slightly more favored than MLR in the day-long problems, and performs signicantly better for the week-long instances. A similar outcome is also observed for the realistic instances, although, for these instances, CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 24 Branch-and-Bound Frac. Vars. (%) jTj Cuts/Nodes Time LP/Root Relaxation OAV 24 4,499 / - 264.1 2.19 / 0.15 168 15,605 / - [0.5%] 1.54 / 0.37 STF 24 478 / 10 53.6 1.09 / 0.13 168 1,940 / - 679.8 0.73 / 0.18 MLR 24 2,914 / - 66.9 2.39 / 0.20 168 11,566 / - 1,243.7 1.76 / 0.29 Table 2.6: Performance measures of branch-and-bound and the % of fractional variables (jGj = 939, Relative optimality gap is reported {in brackets{ when the time limit is hit) the performance delivered by STF is consistently better. Notice that the computational gains observed for the week-long realistic instances are much larger than the total gains observed for the 20 week-long synthetic instances, indicating the potential of STF for instances with more realistic and diverse problem parameters (e.g., distinct generator capacities, ramping rates, or variable-cost functions with many more pieces). Synthetic Instances (summed over 20 instances) Realistic Instances jTj = -5% -2.5% -1% 1% 2.5% 5% -5% -2.5% -1% 1% 2.5% 5% OAV 24 4407.2 4662.3 4363.8 3699.6 4177.9 5343.6 280.3 328.6 233.8 294.5 167.3 254.7 168 1340.2 1154.8 1210.1 1538.2 1484.0 1478.8 [0.5%] [0.7%] [0.5%] [0.8%] [0.4%] [0.4%] STF 24 397.0 367.4 371.8 414.6 491.7 597.5 30.4 33.8 35.8 36.0 35.8 33.9 168 301.2 287.3 300.6 319.8 371.7 348.2 835.7 523.2 761.8 712.5 497.2 658.0 MLR 24 394.3 441.2 367.1 451.1 546 579.3 58.9 47.7 43.5 49.5 37.7 38.4 168 730.4 757.3 1111.0 1293.1 916.2 904.0 1417.7 1847.9 1329.7 1336.2 1362.2 1338.1 Table 2.7: Solution times for modied instances (jGj is 10-fold whenjTj = 24 for synthetic instances, andjGj = 939 for realistic instances). In our nal analysis, we demonstrate the value of considering longer planning horizons in UC problems. Single- day problems are not designed to capture the day to day dynamics of a continuously-evolving demand pattern, and may lead to suboptimal decisions at later periods of the planning horizon. We show this on the realistic instances, by solving 7 sequential day-long UC problems, and their week-long counterpart. Figure 2.2 illustrates the optimal production patterns of two generators under the commitment decisions of these runs. The same pattern appears for most nonbase-load generators. The day-long problems recommend turning these generators o at the end of each day. In contrast, the week-long problem acts with better foresight and recommends less-frequent state transitions, leading to more stable commitment schedules. More importantly, the week-long problem reports 1.8% lower operational costs. 0 12 24 36 48 60 72 84 96 108 120 132 144 156 168 Output (MW) 0 100 200 Time (hr) 0 12 24 36 48 60 72 84 96 108 120 132 144 156 168 Output (MW) 0 100 200 Day-long UCs Week-long UC Figure 2.2: Production patterns of two generators based on the commitment decisions of 7 day-long UC problems, and a week-long problem. CHAPTER 2. A STATE-TRANSITION MIP FORMULATION FOR THE UC PROBLEM 25 2.4. Notes In this study, we developed the state-transition formulation for the UC problem. Our formulation replaces the long-established state variables of the UC formulations with state-transition variables. We compare the performance of the new formulation with two benchmark formulations. We observe that the state-transition formulation performs the best especially for long-horizon instances. The proposed formulation naturally includes valid inequalities that have been used to strengthen alternative formulations. The induced network sub-structure of the formulation allows the eect of these valid inequalities to propagate throughout the planning horizon. Our study also demonstrates the added value of considering longer planning horizons in UC problems. In our deterministic setup, we observed that a week-long planning horizon could lead to cost savings reaching up to 1.8%, compared to the current practice of solving day-long UC problems. The reason is, longer planning horizons allow the UC models to capture demand trends (e.g., across a week), in an eective manner. Furthermore, longer planning horizons could be an important tool to mitigate supply-uncertainty, which has been an emerging issue as ever more wind and solar resources are integrated into power systems. Despite imperfect weather forecasts, the predicted dynamics of renewable supply could be utilized well by longer-horizon UC models, which could lead to decisions with better foresight. Chapter 3 The PH-BAB Framework In recent decades, stochastic mixed-integer programming (SMIP) literature has witnessed a urry of progress. The explicit handling of uncertainty paved the way for better decision-making on a wide scale of real-life problems, and attracted substantial interest from the academic community and practitioners. However, the growing popularity of SMIP also highlighted the need for scalable algorithms that can deliver high-quality decisions in short amounts of time, as SMIP problems are notorious for their computational complexities. Signicant progress has already been made toward addressing certain classes of SMILPs, which, not long ago, baed researchers (see §1.2.2 for examples, and Sen (2005) for a review of methodologies). Despite the advances, the need for an ecient SMIP algorithm without restrictive assumptions { especially in the multi-stage case { is yet to be fullled. In addition, very little is known about the broader class of SMICPs. The purpose of this chapter is to ll these gaps in the literature. In particular, we seek to address the following question: Is it possible to leverage the advances in the stochastic convex programming literature in order to develop a solution framework for multi-stage SMICP with relatively weaker assumptions than similar frameworks in the literature? To answer this question, we will rely on the Progressive Hedging (PH) algorithm of Rockafellar and Wets (1991). The specic contributions and features of this chapter are summarized with the following bullets. We develop a novel branch-and-bound framework for using the PH algorithm to compute optimal solutions to multi-stage SMICPs. We name this framework Progressive-Hedging-based branch-and-bound (PH-BAB). Our study involves a thorough treatment of the convergence of the proposed framework and guarantees optimality under mild assumptions. Within our framework, we implement specialized algorithms that can handle two-stage SMICPs with binary rst-stage and arbitrary second-stage variables. This family of problems have, to date, continued to challenge the stochastic programming community, and often arises in important elds such as telecommunications (e.g., server location problems) and power-systems planning (e.g., UC problems with quadratic production costs). Furthermore, the implemented algorithms serve to bridge the gap between resource-directive and price-directive decomposition paradigms. Our framework is tailored to process, mainly, decomposed convex-relaxations of the original SMIP. Owing to this design, our algorithms are likely to possess better computational performance than its contenders which process deterministic MIP subproblems. In fact, this feature of our framework may yield the highest benet on problems involving more arbitrary convex components (e.g., quadratic objectives), as solving mixed-integer convex problems can be substantially more challenging than solving linear ones. A version of this chapter has appeared in Atakan and Sen (2018). 26 CHAPTER 3. THE PH-BAB FRAMEWORK 27 Among our suggested algorithms, the design of a particular one performs partial optimizations of the decomposed convex-relaxations. As a result, we expect our algorithm to spend minimal eort on regions of the feasible set, bearing no promise. In addition, such a design allows our algorithm to be more robust against suboptimal parameter settings of the PH algorithm. We provide extensive computational experiments on a variety of instances, including a family of extensively- studied linear benchmark instances, and their new quadratic extensions. Our experiments reveal that the assertions of convergence for our specialized algorithms can be observed in short amounts of time for standard test-instances involving up to two thousand scenarios. For convenience of reference, we restate an SMICP for this chapter as follows: min X E[F (X) ] subject to X2N\ S s=1 C s ; x s 2R nr Z nz ;8s = 1:::S: (P) For the denitions of the notations used in the above problem, we refer the reader to §1.2.2. The adopted notation is tailored for a multi-stage stochastic programming setting. It also resembles the notation used in Rockafellar and Wets (1991), easing the transition. We remind the reader of our assumptions, which require the set of scenarios S =f1:::Sg to be nite,C s ;8s = 1:::S; to be closed and convex sets, andf s ;8s = 1:::S; to be convex functions. The organization of the rest of this chapter is as follows. In §3.1, we present a summary of the PH algorithm and discuss its use in the literature as a heuristic method for SMIPs. In §3.2, we develop our framework for using PH to compute optimal solutions to SMICPs in a somewhat general setting. In §3.3, we specialize this framework to the aforementioned class of two-stage SMICPs. The developed algorithms allow our framework, still in its infancy, to compete with its contenders. We provide a computational study in §3.4 which demonstrates the performance of our framework. 3.1. Progressive Hedging For convex stochastic optimization problems, Rockafellar and Wets (1991) have proposed the PH algorithm which, after more than two decades, remains to be one of the popular approaches to multi-stage stochastic opti- mization. In this section, we give a brief overview of their algorithm, along with a discussion of its use for SMIPs. The Algorithm In most stochastic programming problems, the nonanticipativity constraints tie large numbers of scenario subproblems to each other. A natural direction is to dualize these constraints, after which, the subproblems can be independently optimized. However, the coordination of the dual variables may not be trivial, and the optimal dual prices may be dicult to obtain. As a result, the presumed benet of decomposition may never be gained in practice through naive algorithms. Several strategies can be pursued to eciently coordinate the dual variables. For instance, the dual decomposition of Care and Schultz (1999) uses a bundle method (e.g., Helmberg, 2000), whereas Lulli and Sen (2004) adopt a branch-and-price framework. In their PH algorithm, Rockafellar and Wets (1991) address this concern by using a combination of projection and proximal point algorithms (Rockafellar, 1976). We begin the discussion by introducing the multipliers s ;8s = 1:::S, and the associated mapping :s! s . Multiplying the outcomes of this mapping with the scenario probabilities, we obtain p s s . These values should be interpreted as the dual multipliers of the nonanticipativity constraints associated with scenario s. We denote a nonanticipative policy with ^ X : s! ^ x s . Using ^ X, the nonanticipativity restrictions can be explicitly stated as N =fX j X = ^ Xg. Accordingly, the ordinary Lagrangian, achieved through the dualization of the constraint X ^ X = 0, can be dened as follows: L(X; ^ X; ) = S X s=1 p s f s (x s ) + s (x s ^ x s ) , L(X; ) = S X s=1 p s f s (x s ) + s x s ; 82Y: The equivalence of these two denitions stems from the dual feasibility requirements embedded inY, which itself CHAPTER 3. THE PH-BAB FRAMEWORK 28 x 11 =...= x 51 p1 + p2 + p3 p1 p2 p3 p4 p5 p4 + p5 p1 p2 p3 p4 p5 x 12 =x 22 =x 32 t=1 t=2 t=3 t=1 t=2 t=3 t=1 t=2 t=3 t=1 t=2 t=3 t=1 t=2 t=3 t=1 t=2 t=3 P 5 s=1 p s s1 =0 P 3 s=1 p s s2 =0 s=1 s=2 s=3 s=4 s=5 Figure 3.1: The scenario tree representation of a three-stage problem with ve scenarios (left) and an equivalent lifted formulation for scenario decomposition, based on J ornsten et al. (1985) and Guignard and Kim (1987) (each circle denotes a decision). is a consequence of the projected optimality condition 02@L ^ X (X; ^ X; ), i.e., the subgradient ofL ^ X (X; ^ X; ) with respect to ^ X while keeping the outcomes of the maps (X; ) xed. More specically, we have Y = E[ t jA ] = P s2A p s st P s2A p s = 0; 8A2F t ; t = 1:::T ; where, t represents the multipliers corresponding to the nonanticipativity constraints of time period t. Then, for any 2Y and (nonanticipative) ^ X, we have E[ t ^ X t ] =E E[ t ^ X t jA ] =E ^ X t E[ t jA ] = 0: Therefore, when 2Y, the termsp s s ^ x s ;8s, can be omitted from the denition of L(X; ^ X; ). Following this Lagrangian dualization scheme leads to a scenario-wise decomposition of a stochastic program, which we illustrate in Figure 3.1. Let k be the iteration index, and let the superscript k denote the solution of the corresponding variable at the end of the k th iteration. We describe an iteration k of the PH algorithm using the following function: S X s=1 p s f s (x s ) + s (x s ^ x s ) + 2 kx s ^ x s k 2 2 1 2 k s k1 s k 2 2 : (3.1) This function can be interpreted as an extension ofL(X; ^ X; ) with two quadratic regularization terms, each having diering objectives. The termskx s ^ x s k 2 2 ;8s induce better consensus decisions at the end of an optimization, since ^ x s involves components which are shared among certain subsets of scenarios, due to nonanticipativity requirements. The termsk s k1 s k 2 2 ;8s, on the other hand, are proximal penalties which introduce stability to the dual iterates by penalizing their deviations between consecutive iterations. The parameter is a positive constant that scales the impact of these regularizers. An undesired feature of (3.1) is the presence of ^ x s inside the quadratic terms, which destroys the separability of this function with respect to scenarios. The PH algorithm handles this issue by successively projecting two variables and optimizing the third in (3.1). Accordingly, the optimizations with respect to the components of (X; ^ X; ) can be performed independently across the scenarios. Notice that this is the same pattern followed in the broader genre of algorithms, known as the alternating direction method of multipliers (ADMM) (see, for instance, Boyd, 2010). Indeed, the PH algorithm can be seen as a special instance of an ADMM 1 . An iteration k of the PH algorithm begins with xing ( ^ X; ) to their most recent solutions ( ^ X k1 ; k1 ), and solving the following modied subproblems for X: min xs n f s (x s ) + k1 s x s + 2 kx s ^ x k1 s k 2 2 x s 2C s o ; 8s = 1:::S: (P s ) 1 An alternative statement of the problem is given in Rockafellar and Wets (1991). The authors pose the problem as a projected saddle function that is extended with proximal terms on the ^ X and variables. The proof of convergence for the PH algorithm relies on this interpretation. CHAPTER 3. THE PH-BAB FRAMEWORK 29 Above, ^ X k1 can be perceived as a (nonanticipative) incumbent, with respect to which, the variablesX are regular- ized. From a computational perspective, this is generally the most expensive step of the PH algorithm, as it requires the optimization of S convex programs. Next, X is xed to the solution obtained from the modied subproblems X k , is kept at k1 , and the updated form of (3.1) is optimized with respect to ^ X. The optimal solution of this minimization admits a closed form solution, which is given by ^ x k st =E[X k t jA ]; 8s2A; A2F t ; t = 1:::T; where,E [X t jA ] = P s2A p s x k st P s2A p s : (3.2) It is easy to see that the above operation is a projection from the set of arbitrary policies to the set of nonanticipative solutions. Therefore, the resulting policy ^ X k will always be nonanticipative. However, this should not necessarily suggest that ^ X k is feasible with respect toC s ;8s. The nal step of an iteration is the update of the dual multipliers, which is performed by xing (X; ^ X) to (X k ; ^ X k ) and maximizing (3.1) with respect to . Then, we have k s = k1 s +(x k s ^ x k s ); 8s = 1:::S: (3.3) This completes one iteration of the PH algorithm. Observe that individual steps of the algorithm can be performed in parallel, which can be exploited to obtain computationally-powerful implementations. The PH algorithm is reiterated until some termination criteria, such as P S s=1 p s kx k s ^ x k s k 2 <" ph for a suciently small but positive " ph is achieved. It is common practice to set x 0 s as a minimizer in f s (x s ), since by doing so, we also get an immediate lower bound P S s=1 p s f s (x 0 s ) to the optimal objective value of the original problem. Similarly, letting 0 s = 0;8s, leads to 0 2Y, and by induction, one can show that k will then belong toY for any arbitraryk. Finally, usingx 0 s ;8s, one can also determine an initial incumbent ^ X 0 through the procedure in (3.2). A pseudocode of PH is provided in Algorithm 1. Algorithm 1: Progressive Hedging 1 Input: ^ X 0 ; X 0 ; 0 , " ph > 0. 2 Let k = 1 3 while P S s=1 pskx k1 s ^ x k1 s k<" ph do 4 SolvePs for each s = 1:::S, and obtain X k . 5 Form the nonanticipative incumbent ^ X k according to (3.2). 6 Update the dual multipliers according to (3.3). 7 Let k =k + 1 8 end The proof of convergence for the PH algorithm relies on the theory of proximal point algorithms. For details, we refer the reader to Rockafellar and Wets (1991). In what follows, we summarize two important results which will be necessary for the development of our algorithm. We begin with an assumption. A1. All modied subproblems (P s , 8s = 1:::S) encountered during the course of the PH algorithm can be optimized in nitely-many steps. Problems involving linear or quadratic objectives and polyhedral feasible sets clearly satisfy this requirement. Under this assumption, the main convergence result of the PH algorithm is given below. Theorem 3.1 (Rockafellar and Wets (1991)). Consider a convex stochastic optimization problem which satises A1 and the following constraint qualication condition: The only satisfying s 2N Cs (x s ) 8s = 1:::S is = 0; (Q) where N Cs (x s ) is the normal cone associated with the set C s . Assume that the problem has an optimal solution ( ^ X ; ), which is not necessarily unique. Letf ^ X k ; k g 1 k=1 be the sequence generated by Algorithm 1 starting at an CHAPTER 3. THE PH-BAB FRAMEWORK 30 arbitrary point ( ^ X 0 ; 0 ). Then, when k!1, we have ^ X k ! ^ X and k ! : Furthermore, in terms of the norm expressionk(X; )k r = (kXk 2 + 2 kk 2 ) 1 =2 , one will have in every iteration k that k( ^ X k+1 ; k+1 ) ( ^ X k ; k )k r k( ^ X k ; k ) ( ^ X k1 ; k1 )k r : (3.4) Proof. See Theorem 5.1 of Rockafellar and Wets (1991). Remark 3.1. The constraint qualication condition (Q) is super uous when the considered problem involves linear or quadratic objective and polyhedral feasible set (see Theorem 4.3 in Rockafellar and Wets (1991)). The existence of nite algorithms for arbitrary families of convex optimization problems is not evident. In such cases, Rockafellar and Wets (1991) suggest optimizingP s in an approximate sense. Following their scheme, A1 can be relaxed and, under fairly general conditions, the convergence result in Theorem 3.1 will remain valid (see Theorem 5.4 and the related discussion in Rockafellar and Wets (1991)). For this reason and to keep this paper concise, A1 will be presumed throughout our study. Notice that an optimal solution pair ( ^ X ; ) leads to an optimal solution X , since the condition X = ^ X is necessary for the PH algorithm to terminate. The following theorem provides a saddle point characterization for an optimal solution (X ; ), obtained by the PH algorithm, to the underlying stochastic convex optimization problem. Theorem 3.2 (Rockafellar and Wets (1991)). The optimal solution pair (X ; )2NY obtained by the PH algorithm is a saddle point of the ordinary Lagrangian L(X; ), relative to minimizing over x s 2C s ;8s = 1:::S, and maximizing over 2Y. Proof. See Theorem 4.2 of Rockafellar and Wets (1991). Use of PH on SMIPs Contrary to the case of convex stochastic optimization problems, strong theoretical guarantees on the convergence of the PH algorithm for non-convex problems are not available. In the case of SMIPs, the modied subproblems solved at each iteration become non-convex because the feasible set ofP s (C s \(R nr Z nz )) is discrete. Empirical evidence indicates that executing the PH algorithm with MIP subproblems may not always lead to an optimal solution, and can even lead to oscillatory behavior of the iterates (Gade et al., 2016). However, such diculties can be overcome with algorithmic novelties, and the PH algorithm can be used as a viable heuristic approach for solving SMIPs. For a variety of SMILPs, it has been shown that the PH algorithm indeed provides promising solutions (see, for instance, Takriti et al., 1996; Lkketangen and Woodru, 1996; Ryan et al., 2013). Nevertheless, due to non-convexity, the convergence of the algorithm to a (near)-optimal solution relies heavily on the quality of the implementations rather than mathematical assertions. In the remainder of this paper, we will refer to similar uses of the PH algorithm on SMIPs as PH-Heuristic. Several studies aim to improve the computational eectiveness of PH-Heuristic. For instance, Watson and Woodru (2010) analyze the eect of the penalty parameter on the speed and the solution quality of the algorithm. Another line of research aims to reduce the duality gap that arise due to the relaxation of the nonanticipativity constraints in a mixed-integer program. Gade et al. (2016); Ryan et al. (2016) show that enforcing the nonan- ticipativity restrictions for subsets of scenarios (scenario bundling) is shown to improve the quality of solutions. While not directly related to PH-Heuristic, studies on eliminating the duality gap through augmented Lagrangians (Boland and Eberhard, 2014; Feizollahi et al., 2017) may carry signicant potential in providing a theoretical basis for such algorithms. For the purpose of achieving convergence guarantees, Guo et al. (2015) hybridize the PH-Heuristic with the dual decomposition of Care and Schultz (1999). Boland et al. (2016) propose an extension of the PH algorithm, which CHAPTER 3. THE PH-BAB FRAMEWORK 31 converges to the optimal objective value of the Lagrangian dual problem, achieved through the dualization of the nonanticipativity constraints inP. Their study incorporates a simplicial decomposition method (SDM) into the PH algorithm. The SDM successively provides inner-approximations of the convex hull of C s \ (R nr Z nz ), therefore resolves the non-convexity and leads to a convergent algorithm. The suggested method is designed to solve one MIP subproblem per scenario, therefore achieves somewhat similar computational performance as the aforementioned uses of the PH algorithm. Finally, the recent work of Barnett et al. (2017) propose to use a branch-and-bound procedure as a wrapper for the PH-Heuristic. Their framework uses the incumbents ( ^ X k ) to partition the feasible region over which PH-Heuristic is iterated. The components of this partition are fathomed based on Lagrangian bounds on their optimal objective values, which are obtained by solving S mixed-integer subproblems (see Gade et al., 2016). The provided computational experiments reveal that the suggested framework spends substantially more computational eort than the standard PH-Heuristic in order to achieve higher-quality solutions, however lags somewhat behind the performance of a commercial mixed-integer programming solver that is directly applied toP. These outcomes might be interpreted as symptoms of heavy-reliance on mixed-integer programming technology in iterative decomposition algorithms. We emphasize that the PH-Heuristic algorithm and its aforementioned extensions require the solution of S mixed-integer programming subproblems at each iteration. In general, a stochastic program is only meaningful when a large-enough set of scenarios are considered. Depending on how large the sample size is, the computational eort expended on iterating the PH algorithm may actually defeat the purpose of decomposition. We intend our framework to remedy this. 3.2. General Framework of PH with Branch and Bound The power of branch-and-bound strategies in dealing with the non-convex feasible sets of MIPs is undisputed. In our framework, we integrate their potential with the scenario-decomposition capability of the PH algorithm. The resulting framework, which we refer to as PH-BAB, can be perceived as a branch-and-bound method that uses the PH algorithm to solve the nodal relaxations, approximately (see Figure 3.2). Start with the convex-relaxation of SMICP Update the incumbent & perform fathoming Use PH to solve the nodal relaxation (initialize with parent’s ˆ X,⇤ ) Branch on a fractional variable false true Nodes left? integer? true Return Solution false Select a branch-and-bound node Figure 3.2: A basic PH-BAB algorithm. Let Q s =fx s 2R n j l s x s u s g denote the feasible set dened by bounds on the vector x s . The mapping Q :s!Q s will be convenient to collectively represent the bounds across all scenarios. Accordingly, at any given branch-and-bound node , the corresponding mapping will be represented with Q : s! Q s . The feasible sets of individual scenario subproblems at node , i.e., C s \Q s , will be denoted with the short-hand notation C s . CHAPTER 3. THE PH-BAB FRAMEWORK 32 Whenever we say that node is being solved, we mean that the following nodal relaxation is being optimized: z() = min X E[F (X) ] subject to X2N\ ( S s=1 C s ): (R ) The root node of the branch-and-bound tree, denoted witho, will possess the bound constraintsQ o . Without loss of generality, we assume that C s C o s ;8s = 1:::S. Therefore, nodeo should be interpreted as the convex relaxation ofP, achieved by removing the integrality restrictions. The pseudocode of a basic PH-BAB algorithm is provided in Algorithm 2. Algorithm 2: A Basic PH-BAB Algorithm 1 Initialization: Create the root node o, letH fog. Set zU = +1, zL =1, ^ XU = NULL. 2 whileH6=; do 3 Pick node fromH, and letH Hnfg. 4 Determine the feasibility of . 5 if is feasible then // else, is fathomed 6 Invoke the PH algorithm to solve . 7 Get solution ^ X of , and its optimal objective value z(). 8 if z()<zU then // else, is fathomed 9 if ^ X is fractional then 10 Pick a fractional component ^ x stj of ^ X . 11 Create nodes + and . Let Q + =Q =Q , while introducing the additional restrictions lstj =d^ x stj e;8s2FFt to Q + , and ustj =b^ x stj c;8s2FFt to Q . 12 LetH H[f ; + g. 13 else 14 Let ^ XU = ^ X , and zU =z(). 15 Fathom:H Hnf 0 jz( 0 )zUg. 16 end 17 end 18 end 19 Let zL = minfz()j2Hg. 20 end We draw important distinctions on the use of the branch-and-bound method among the PH-BAB and PH- Heuristic algorithms. Both algorithms leverage the capability of branch-and-bound for handling non-convex feasible sets, however, their goals are dierent. In PH-BAB, we seek an optimal solution toP, and form a single branch- and-bound tree for the feasible set of the entire problem. In the latter, branch-and-bound trees are formed to nd optimal solutions to modied subproblems, which are modeled as SMICPs. In fact, in the PH-Heuristic algorithm, the branching decisions are performed on X variables { independently { across the scenarios. As a result, the optimal solutions for the modied subproblems need not agree with each other with respect to the nonanticipativity requirements, and may even produce incumbents ^ X that disobey integrality restrictions (e.g., consider the averaging of binary variables). Moreover, after an iteration is completed, the branch-and-bound trees are usually discarded, and information are not conveyed from one iteration to another. In contrast, the single branch-and-bound tree of the PH-BAB algorithm progressively partitions the feasible set ofP into ner feasible regions. The branching essentially occurs on ^ X, which is re ected on X by simultaneously branching on the associated observationally- indistinguishable decisions (see line 11 of Algorithm 2). With this branching strategy, solutions that violate the nonanticipativity restrictions are progressively discarded from future consideration. As a result, as the depth of the branch-and-bound tree grows, the nonanticipativity requirements imposed on the integer variables can be attained with greater ease. Ecient implementations of PH-BAB easily exploit this feature by giving branching priorities on the variables in the increasing order of their time periods. As in any other branch-and-bound-based method, the PH-BAB can signicantly benet from warm starts when solving its nodal relaxations. Indeed, the initialization of the PH algorithm that is suggested in §3.1, should only be invoked at the root node, as it involves the optimization of S minimization problems. For all other branch-and- CHAPTER 3. THE PH-BAB FRAMEWORK 33 bound nodes, the nal primal and dual solutions of the parent node can be used as a starting point. For faster primal convergence, the primal solution inherited from the parent node should be projected onto the set dened by the new variable bounds. These simple considerations can carry signicant potential in improving the performance of PH-BAB implementations. We leave further computational details to later sections, and continue with a discussion on the convergence of the proposed method. The following assumptions will be necessary for this purpose: A2. All integer variables in problemP have nite lower and upper bounds. A3. Each feasible nodal relaxation satises the constraint qualication conditionQ. A4. A certicate of feasibility or infeasibility for each nodal relaxation (i.e., for the setN\ S s=1 C s ) can be computed in nite time. A5. For each feasible node, the PH algorithm yields an optimal solution in nite time. The presence of A2 does not restrict the practicality of our framework as arbitrarily loose bounds can be enforced on unbounded integer variables. Assumption A3 is a technical condition that is satised in most problem formulations, and holds when the setsC s ;8s, are polyhedral (see Remark 3.1). Assumption A4 may lead to concern as infeasibility detection for arbitrary convex problems is not trivial. For well-understood families of convex optimization problems (e.g., linear or quadratic programs), the state-of-the-art solvers are equipped with mechanisms which can detect the consistency of C s and determine the feasibility of the corresponding subproblem. Pathological cases may still occur when C s 6=;;8s = 1:::S, butN\ S s=1 C s =;, in which case the infeasibility must be detected with alternative means, for instance, by considering a feasibility problem before decomposing the nodal relaxation. For most practical purposes, such issues seldom occur, and therefore will not be elaborated further, but instead, left for future investigation. While these three assumptions will be necessary throughout this paper, we state A5 exclusively for the next proposition. Proposition 3.1. For problemP which satises A2-A5, Algorithm 2 yields an optimal solution in nite number of steps. Proof. Algorithm 2 diers from a standard branch-and-bound algorithm only with respect to its simultaneous branching rule on observationally-indistinguishable decisions. However, in a feasible solution, these variables must assume equal values due to nonanticipativity restrictions, therefore the branching rule can only discard nodes that are surely infeasible. Accordingly, the proof of convergence follows that of a standard branch-and-bound algorithm, performing exact evaluations. For a rigorous proof, we refer the reader to Horst and Tuy (2013), but we provide the basic argument below. At every iteration of the algorithm, a single node is removed fromH for processing. By A4, an infeasible node can be detected and discarded in nite time. Due to A5, a feasible node can be optimized by the PH algorithm in nite time. If this optimization leads to an integer-feasible solution, the cardinality ofH must reduce in the next iteration, asH will no longer contain . If the optimization produces a fractional solution, branching will occur andH will gain one more element. Due to A2, the feasible set of the problem can only contain a nite number of integer points. As a result, branching can only occur nitely-many times, thereforeH must empty in a nite number of iterations. At this point, the incumbent solution will be optimal toP. A critical concern can be raised to A5, as the nite termination of the PH algorithm is not guaranteed by the assertion given in Theorem 3.1. Alternative PH-BAB implementations can be studied to ensure global convergence. In the rest of this section, we study a particular design where the nodes of the branch-and-bound tree are partially solved, and fathoming is performed according to bounds on their optimal objective values. In particular, after each (partial) optimization of node , a lower bound (z L ()) and an upper bound (z U ()) will be computed on its optimal objective value (z()). In our setup, these bounds will be derived from the solutions generated by the PH algorithm. We require that, when necessary, these bounds can be made arbitrarily tight, so that z() can be determined by observing the equality ofz L () andz U (), perhaps, in the limit. We begin with our lower-bounding CHAPTER 3. THE PH-BAB FRAMEWORK 34 procedure, which utilizes the solution for at an arbitrary iteration of the PH algorithm. The following proposition is based on Lagrangian duality, and its proof follows Gade et al. (2016). Proposition 3.2. Consider a feasible nodal relaxationR with an optimal solution ^ X , and assume it satises A1. Suppose that the PH algorithm is used to solve the problem. Let the multipliers k , obtained at the end of the k th iteration of the algorithm, satisfy k 2Y. Dene the function D s ( k s ) as D s ( k s ) = min xs f s (x s ) + k s x s jx s 2C s 8s = 1:::S: (3.5) Then,D( k ) = P S s=1 p s D s ( k s ) is a lower bound on the optimal objective function valuez() of the nodal relaxation. Proof. Since each subproblem in (3.5) is minimized with respect to x s , the following inequality is true: D( k ) S X s=1 p s f s (^ x s ) + k s ^ x s = S X s=1 p s f s (^ x s ) + S X s=1 p s k s ^ x s =z() +E[ k ^ X ]: (3.6) Furthermore, for any time period t, we have E[ k t ^ X t ] =E E[ k t ^ X t jA ] =E ^ X t E[ k t jA ] (since ^ X 2N ) = 0: (since k 2Y) Applying this argument for every t leads the expectation in (3.6) to vanish, which completes the proof. It should be noted that the Lagrangian bound suggested in Gade et al. (2016) and adopted by Barnett et al. (2017) replaces the convex set C s with C s \ (R nr Z nz ) in (3.5). As mentioned earlier, our framework intends to avoid non-convexity in its subproblems, even though the bounds that they provide may be signicantly better. A6. For any feasible node with an optimal objective value z(), there exists a nite value U z(). The above assumption is not essential for the convergence of our algorithm, however it ensures that the fathoming in the branch-and-bound tree is not delayed for innitely-many steps on any node. Our upper-bounding procedure, which utilizes the solutions to ^ X at any iteration of the PH algorithm, is provided below. Proposition 3.3. Consider a feasible nodal relaxationR . Assume that it satises A1, and admits an upper bound U , as dened in A6. Suppose the PH algorithm is used to solve this problem. Let ^ X k denote the incumbent obtained at the end of the k th iteration of the algorithm. Dene V s (^ x k s ) as follows: V s (^ x k s ) = ( f s (^ x k s ) for ^ x k s 2C s U otherwise. Then, V ( ^ X k ) = P S s=1 p s V s (x k s ) is an upper bound on the optimal objective value z() of the nodal relaxation. Proof. Since V ( ^ X k ) simply evaluates the objective value of the (nonanticipative) policy ^ X k , the statement is obviously true. Remark 3.2. The a priori knowledge of U for all may appear restrictive, however, within a branch-and-bound procedure, they can be replaced with a nite global upper bound on the optimal objective value of problemP (e.g., through a feasible-solution heuristic). This may lead to situations where V ( ^ X k )<z() for some node , however, such nodes can never contain a better integer-feasible solution than the one on hand, and should be fathomed in any case. CHAPTER 3. THE PH-BAB FRAMEWORK 35 The following theorem formalizes the convergence condition on the suggested bounds that is sought in our analysis. Theorem 3.3. Consider a feasible nodal relaxationR . Assume that it satisesA1 and the constraint qualication conditionQ, and admits an upper bound U , as dened in A6. Let ( ^ X ; ) denote an optimal solution to this problem, and z() represent its optimal objective value. Suppose the PH algorithm is used to solve this problem. Consider the sequence of lower bounds D( k ) 1 k=1 and upper bounds V ( ^ X k ) 1 k=1 generated at the end of the k th iteration of the PH algorithm, according to Propositions 3.2 and 3.3, respectively. Then, for k!1, we have D( k )z()V ( ^ X k ); and D( k )!z(); V ( ^ X k )!z(): (3.7) Proof. LetL(X; ) represent the ordinary Lagrangian for node , and deneI as 0 if ^ x s 2C s ;8s = 1:::S, and U if otherwise. ( ^ X) = ( 0 ^ x s 2C s \Q s ;8s = 1:::S; U otherwise. Then, for any iteration k, the assertions of Propositions 3.2 and 3.3 can be stated as follows: D( k ) = min xs2C s ;8s L(X; k )L( ^ X ; k ) L( ^ X ; ) =z() =L( ^ X ; ) +I ( ^ X ) since ^ x s 2C s ;8s L( ^ X k ; k ) +I ( ^ X k ) since ^ X2N andI ( ^ X k ) 0 =V ( ^ X k ): This proves the inequality in (3.7). By Theorem 3.1, we have ( ^ X k ; k )! ( ^ X ; ) as k!1, implying that L( ^ X k ; k )!L( ^ X ; ). By Theorem 3.2, ( ^ X ; ) is a saddle point of L( ^ X; ), therefore ^ X 2N\ S s=1 C s , and L( ^ X ; ) = min xs2C s ;8s L(X; ) =D( ) =V ( ^ X ). As a result, as k!1, we have V ( ^ X k )D( k )! 0, which completes the proof. Propositions 3.2 and 3.3 can be invoked at any arbitrary iteration of the PH algorithm, even before it approaches to an acceptable solution. The convergence of the PH algorithm can be monitored by periodically evaluating these bounds, and a solution of a node can be labeled as optimal when the gap between these bounds fall within numerical tolerances. When a node is only partially solved, fathoming of the branch-and-bound tree can be performed according to these lower bounds. If a (partial) optimization produces a primal solution that satisfy the integrality requirements, then the upper bounding procedure will determine the objective value of a candidate solution to the original problem (despite the fact that the node may still need to be solved to optimality). Alternative bounding procedures can also be investigated. For instance, in Proposition 3.2, f s can be replaced with its minorants, and in Proposition 3.3, special problem structures can be exploited. However, these alternatives must satisfy the convergence condition (3.7), stated in Theorem 3.3. The following specications must be fullled in the design of our next algorithm. A7. The branch-and-bound node selection is performed according to the least lower bound rule. A8. When the PH algorithm (partially) solves a node for the i th time, it is terminated, either, (a) according to a nite iteration limit, or (b) with a convergence tolerance " i ph > 0 where " i ph ! 0 as i!1, or (c) with an optimality tolerance i ph > 0 where i ph ! 0 as i!1. Furthermore, the next optimization of is initialized with the nal solutions of (X; ^ X; ), obtained at the end of the i th run (along with all other possible algorithmic settings of this run). A9. For each node , the computed lower bounds (z L ()) and upper bounds (z U ()) satisfy the convergence condition (3.7) given in Theorem 3.3. CHAPTER 3. THE PH-BAB FRAMEWORK 36 The latter two assumptions ensure that the PH algorithm can continue solving each nodal relaxation seamlessly, despite being interrupted after a nite number of steps. As a result, at each node, the convergence of the PH algorithm to an optimal solution will follow Theorem 3.1, and the convergence of the bounds z L () and z U () will follow Theorem 3.3. We present a PH-BAB algorithm that performs partial solves of its nodes in Algorithm 3. In the ensuing proposition, we show that implementing Algorithm 3 according to A7-A9 leads to an exact method for optimizing multi-stage SMICPs. Algorithm 3: A PH-BAB with Partial Node Solves 1 Initialization: Create the root node o, letH fog. Set zU = +1, zL =1, ^ XU = NULL. 2 whileH6=; do 3 Pick node fromH and letH Hnfg. 4 Determine the feasibility of the sets C s ;8s = 1:::S. 5 if C s ;8s; is feasible then // else, is fathomed 6 Invoke PH to (partially) solve . 7 Get solution ^ X of , and the lower bound zL(). 8 if zL()<zU then // else, is fathomed 9 if ^ X is fractional then 10 Pick a fractional component ^ x stj of ^ X . 11 Create nodes + and . Let Q + =Q =Q , while introducing the additional restrictions lstj =d^ x stj e;8s2FFt to Q + , and ustj =b^ x stj c;8s2FFt to Q . 12 LetH H[f ; + g. 13 else 14 Obtain an upper bound zU () and the corresponding solution ^ X U . 15 if zU ()<zU then 16 Let ^ XU = ^ X U , zU =zU (). 17 Fathom:H Hnf 0 jzL( 0 )zUg 18 end 19 LetH H[fg. // is replaced 20 end 21 end 22 end 23 Let zL = minfzL()j2Hg. 24 end Proposition 3.4. Consider problemP which satises A1-A4, and let Algorithm 3 be implemented according to A7-A9. Assume that the problem has at least one optimal solution with the objective value z ? . Dene j as the iteration index of Algorithm 3, and let z j L 1 j=0 and z j U 1 j=0 denote the sequences of best lower and upper bounds produced according to lines 23 and 16, respectively. Also, let X j U 1 j=0 be the sequence of candidate solutions (X U ) produced by the algorithm in line 16. Then, as j!1, we have z j L %z ? and z j U &z ? : (3.8) Furthermore, every accumulation point of the sequence X j U 1 j=0 will be an optimal solution toP. Proof. Due to A4, all infeasible nodes that are encountered by the algorithm will be discarded in nite time, and due to A8, all feasible nodes that are encountered will be (partially) optimized in nite time. As a result, Algorithm 3 expends only nite amount of time on each node that it decides to process. In fact, by design, the algorithm performs branching after each (partial) optimization of a feasible node, when the node produces a fractional solution. Due to A2,9 a nite iteration after which the branching process must stop producing nodes. From this iteration forward, all nodes that stay inH must satisfy the integrality restrictions. Without loss of generality, the rest of the proof will only consider the nodes that are processed subsequently. For the rest of the proof, we deneJ to be the set of iteration indices, at which, a particular node is selected and processed by the algorithm. CHAPTER 3. THE PH-BAB FRAMEWORK 37 We begin with showing that z j L 1 j=0 converges to z ? . By way of contradiction, suppose this is not true. This implies that9 a node 2H which determines z j L 1 j=0 innitely-many times, and has the subsequence of lower bounds z j L () j2J not converging to z ? . Due to A9, this subsequence must have the limit z(). Suppose z() =z ? +c. Clearly,c< 0 contradicts the optimality ofz ? , leaving the only possibilityc 0. Forc> 0, we must have z j 0 L () z ? +" for nite j 0 and "2 (0;c), since the subsequence z j L () j2J can only have nitely-many elements in any positive neighborhood of z ? +c. However, at least one other node ? 2H must contain an optimal solution, and for this node, we can write z j 0 L ( ? ) z ? < z ? +" z j 0 L (). As a result, A7 ensures that node cannot determine z j L when jj 0 , unless c = 0. In fact, if c> 0, an upper-bounding sequence z j U 1 j=0 !z ? will ensure that such will be removed fromH after a nite number of iterations. We conclude that, only the nodes with the optimal objective value z ? can determine the best lower bound for innitely-many times, and for such , A9 ensures that lim j!1;j2J z j L () = z ? . Since z j L = z j L (),8j2J and 2H, we have every subsequence of z j L 1 j=0 converging to z ? , therefore lim j!1 z j L =z ? . Next, we show the convergence of z j U 1 j=0 to z ? . Since all 2H satisfy the integrality restrictions, line 14 of Algorithm 3 will be invoked at every iteration. Moreover, due to A7 and the above discussion, all nodes that the algorithm processes innitely-many times, will have the optimal objective value z ? . For any such , A9 ensures that lim j!1;j2J z j U () = lim j!1;j2J z j L () =z ? . Since z j U = min z j U ()j2H , we can only have lim j!1 z j U = z ? without invalidating the optimality of z ? . As a result of this (and the previous discussion), the assertion in (3.8) must be true. Finally, our claim on the sequence X j U 1 j=0 can be trivially conrmed by noting that the best upper boundz j U , produced at iterationj, correspond to the objective value of the feasible solution ^ X j U , which, by virtue of Theorem 3.1, converges to some optimal solution ofP. Proposition 3.4 suggests that Algorithm 3 can resolve the complexity caused by the discrete nature of the feasible set ofP, in nite time, and thereafter, convergence to an optimum must happen. Observe, also, that assumption A6 plays an implicit role in the course of the algorithm, by guaranteeing that all suboptimal nodes will be discarded after nitely-many steps. Remark 3.3. Algorithm 3 always performs branching on nodes whose (partial) optimizations produce fractional solutions. Depending on the application, it might be benecial to delay branching until the node at hand is suciently explored. In such designs, the assertion in Propositions 3.4 will remain valid, provided that branching on any node is never postponed for innitely-long. 3.3. PH-BAB for a Special Class of Two-Stage Stochastic Programs Similar to all branch-and-bound-based methods, the PH-BAB framework relies on enumeration, and straight- forward implementations are unlikely to provide promising results for large-scale problems. Computational per- formance can be enhanced through the use of cutting planes, improved branching and node-selection rules, and many other adjunct procedures, all of which have proved themselves in the eld of deterministic mixed-integer pro- gramming. Moreover, it is well known that special problem structures, which commonly appear in many stochastic programming formulations, can lead to powerful algorithms. The goal of this section is to make use of these facts and develop a computationally eective PH-BAB implementation. The resulting implementation will resemble the branch-and-x coordination of Escudero et al. (2009) (see, also, Pages-Bernausa et al., 2015). We will focus on a special problem structure, which allows us to directly utilize the capabilities of state-of-the-art mixed-integer programming solvers for eciently solving MICP subproblems. The structure we consider is a two-stage stochastic program with only binary variables in its rst-stage. Remarks will be provided for future extensions. For clarity of exposition, we express some of our notation within the two-stage context. We denote the number of continuous and integer variables at t th stage of the stochastic program as n t r and n t z , respectively, and let n t = n t r +n t z . The nonanticipativity restrictions are not meaningful for the second-stage variables and will be CHAPTER 3. THE PH-BAB FRAMEWORK 38 omitted throughout this section. Accordingly, a nonanticipative policy can be constructed from the expectation of the rst-stage decisions x s1 ;8s, alone. We refer to the outcome of this expectation as x. With respect to our previous notation, we have ^ X :s! ^ x s = ( x; x s2 );8s, suggesting that x; ^ x 11 ::: ^ x S1 are all synonyms of each other. The set of nonanticipative solutions can then be represented as follows: N = n (x 11 :::x S1 )2 S s=1 R n 1 jx s1 = x;8s = 1:::S o : This state-vector representation of nonanticipativity allows separable Lagrangians as in Higle and Sen (2006) (see, also Lubin et al., 2013). At each scenario, we partition the feasible region of the problem as C s = (C 1 R n 2 )\C s2 ; where C 1 R n 1 and C s2 R n 1 +n 2 are closed and convex sets. The dropped scenario indices indicate that the corresponding set is not aected by the realization of the randomness. We dene the objective function of the problem as f s (x s ) =f 1 (x s1 ) + min xs2 n f s2 (x s2 )j (x s1 ;x s2 )2C s2 ; (x s1 ;x s2 )2R n 1 R n 2 r Z n 2 z o ; (3.9) where f 1 :R n 1 !R and f s2 :R n 2 !R are both convex functions. The inner minimization in (3.9) is assumed to have the objective value +1 when its input x s1 does not lead to a feasible x s2 . This minimization has no eect on the way the PH iterations are performed, however gives an alternative interpretation on the objective of the problem. Minimizingf s can be perceived as searching for the best rst-stage solution, since the inner minimization in (3.9) (a.k.a., the second-stage subproblem) will be shaped according to its value. Moreover, the altered form off s aects the way the upper bound evaluations in Proposition 3.3 are performed. For instance, for a (nonanticipative) rst-stage decision x produced from some node , the evaluations take the form f 1 ( x ) + S X s=1 p s min xs2 f s2 j ( x ;x s2 )2C s2 ; where C s2 = C s2 \Q s . Accordingly, the procedure essentially assesses the quality of the rst-stage decisions. While computing an upper bound now requires the solution of S second-stage convex subproblems, the likelihood of nding a feasible upper bound is much higher. Following this discussion, we letB =f0; 1g, and pose a two-stage SMICP with pure binary rst-stage variables as follows: min X E[F (X) ] subject to: x s1 = x;8s = 1:::S; x s 2C s ;8s = 1:::S; and x s1 2B n 1 ;8s = 1:::S: (P-2) An alternative representation of (P-2) within the stagewise decomposition context will be useful. In this setting, the rst-stage problem will be unique, and dened as min f 1 ( x) +E[H( x) ]j x2C 1 \B n 1 : (3.10a) Observe that the rst-stage problem is optimized over the (nonanticipative) x variables. The second-stage subprob- lems (i.e.,H( x) :s!h s ( x)) are parameterized according to the solution of x in the above problem, and dened as follows: h s ( x) = min xs2 f s2 (x s2 )j ( x;x s )2C s2 ; x s2 2R n 2 n 2 z Z n 2 z ; 8s = 1:::S: (3.10b) As suggested by the previous discussion, a key characteristic of stochastic programming problems is that for a given nonanticipative solution for the rst-stage, the second-stage subproblems can be evaluated independently. Indeed, stochastic programming algorithms that are in the stagewise-decomposition context mainly rely on this fact. The main purpose of this section will be to investigate under which circumstances the PH-BAB algorithm can also leverage this property, so that promising integer-feasible solutions can be easily produced using the state-of-the-art mixed-integer programming solvers. Despite performing scenario-wise decomposition, the PH-BAB algorithm can eventually end up with nodes, CHAPTER 3. THE PH-BAB FRAMEWORK 39 where the nonanticipativity is guaranteed to hold. Such situations occur when the branch-and-bound procedure restricts the lower and upper bounds of all rst-stage variables to unique points. In these cases, the optimization of a nodal relaxation reduces to solving min xs2 f 1 (x s1 ) +f s2 (x s2 )j (x s1 ;x s2 )2C s2 ; 8s = 1:::S; where, for a xed rst-stage solution x , Q involves the nonanticipativity restriction x s1 = x ;8s. The above subproblems can be independently optimized using any convex programming solver. When f s2 is linear and C s2 is polyhedral, these subproblems turn into linear programs, which can be eectively handled in nite-time. Moreover, from this point forward, the role of the PH algorithm will essentially be over for node and its children. As a result, deterministic mixed-integer programming solvers may now be employed to fulll the integrality restrictions in the second-stage subproblems. Although the benet of a xed rst-stage solution is apparent, whether such situations can be eciently induced must be claried. For problems involving continuous rst-stage variables, our algorithm, in its current form, deems this impossible, and must be extended to allow branching on continuous variables. For pure integer rst- stage variables, we require the branch-and-bound tree to be exhausted for all promising branches, therefore the performance of the resulting algorithm will be unsatisfactory. However, when the rst-stage only involves pure binaries (as we have assumed in this section), an ecient branching scheme can be developed. Our scheme does not disrupt the overall ow of the branch-and-bound procedure, but is invoked only when a node produces an integral rst-stage solution. Consider a node, whose rst-stage variables are not completely xed by the branch-and-bound algorithm. For the sake of argument, suppose that the PH algorithm terminates with an integral rst-stage solution ( x ) but the second-stage is still fractional. For these situations, our recommendation is to follow a modied branching strategy. We create two nodes f and : in such a way that x is not feasible in : , but it is the only feasible rst-stage solution in f . In other words, these two nodes are created according to the unconventional disjunction x s1 = x or x s1 6= x . The new bounds for f can be expressed eortlessly as Q f s =Q s \ n x s 2R n 1 +n 2 j x x s1 x o ; 8s = 1:::S: (3.11) On the other hand, describingQ : s is not straightforward since the imposition x6= x is non-convex. At this point, the pure-binary-variables assumption at the rst stage allows the use of the inequality X j: x j =1 x s1j X j: x j =0 x s1j jfj : x j = 1gj 1; (3.12) which cuts o x from the corresponding feasible set. The above inequality has been commonly used for similar purposes in the stochastic-programming literature (see, for instance, Sen and Sherali, 2006; Ahmed, 2013). With a little stretch on the meaning of Q : ! , we can dene the variable bounds in : as follows: Q : s =Q s \ n x s 2R n 1 +n 2 j (3.12) holds o ; 8s = 1:::S: (3.13) This concludes our branching procedure. As a consequence of the above suggestions, our algorithm can now produce rst-stage solutions, with respect to which, the second-stage value functions can be evaluated. An exact evaluation would require imposing the integrality restrictions on the corresponding second-stage variables, and solving the resulting subproblems using deterministic MICP solvers. In essence, such an operation forms a branch-and-bound sub-tree that emanates from the node involving xed rst-stage solutions. Determining an optimal solution to this sub-tree (if such a solution exists) will produce a candidate optimal solution to the original problem. However, it should be noted that the exact evaluations of the second-stage value functions may be computationally expensive, despite the power of commercial solvers. Therefore, these evaluations should only be pursued when the corresponding rst-stage CHAPTER 3. THE PH-BAB FRAMEWORK 40 solutions are deemed promising. Indeed, instead of frequently invoking MICP solvers, algorithms that successively provide tighter approximation of the second-stage value functions (such as Sen and Higle, 2005; Sen and Sherali, 2006; Gade et al., 2014; Zhang and K u c ukyavuz, 2014) could be preferable. Their potential within the PH-BAB scheme should be investigated in the future. Along these lines, we consider all nodes, which involves xed rst-stage variables x , as MICPs of the following form: f 1 ( x ) + min xs2 n f s2 (x s2 )j ( x ;x s2 )2C s2 ; ( x ;x s2 )2R n 1 R n 2 r Z n 2 z o ; 8s = 1:::S: (3.14) Using a presumably fast MICP solver, we optimize the above subproblems in two phases. First, we perform a partial optimization to obtain lower bounds on the optimal objective value of the subproblems. These lower bounds immediately lead to a lower bound on the objective value of the best integer-feasible solution that can contain (if it contains any). We refer to these bounds as z L (), although this notation should not be confused with the lower bounds generated for other arbitrary nodes according to Proposition 3.2. Using z L (), it is possible to assess whether needs to be fathomed, by comparing it with the best upper bound z U of the PH-BAB algorithm. If the node may indeed contain a better integer-feasible solution, we continue with the second phase of our node algorithm where the subproblems (3.14) are completely optimized. An outline of our procedure is given in Algorithm 4. Algorithm 4: A Node Algorithm 1 if node does not have xed rst-stage variables then 2 Determine the feasibility of . 3 if is feasible then Invoke PH to (partially) solve . // else, will be fathomed 4 else 5 Enforce the integrality restrictions on the second-stage variables to form (3.14). 6 Determine the feasibility of (3.14). 7 if is feasible then // else, will be fathomed 8 Partially solve (3.14), and obtain a lower bound zL(). 9 if zL()>zU then // else, will be fathomed 10 Completely solve (3.14), and obtain the optimal objective value z(). 11 end 12 end 13 end The strategies presented in this section can be incorporated to both Algorithm 2 and 3, in order to increase their computational eciencies. In Algorithm 5, we present the pseudocode for the extension of the latter. In comparison to Algorithm 3, this version of the PH-BAB method can converge to an optimal solution in nite time, provided that the following assumption holds. A10. For any set of problems of the form (3.14), there exists an MICP solver which can either determine the infeasibility of any one of the subproblems, or yields an optimal solution to all, both in nite number of steps. The assumption clearly holds when the feasible sets C s ;8s, are polyhedral, and the functions f s2 ;8s, are linear or quadratic. As a consequence of A10, lines 5-12 of Algorithm 4 can be executed in nite time. Proposition 3.5. Consider problemP-2 which satises A1-A4, A10, and let Algorithm 5 be implemented accord- ing to A7-A9. Then, Algorithm 5 yields an optimal solution in nite number of steps. Proof. The proof is mainly based on the convergence of Algorithm 3, as stated in Proposition 3.4. Indeed, Algorithms 3 and 5 dier mainly with respect to the introduced branching strategy and the node algorithm. Algorithm 5 iterates until all of the rst-stage variables are branched on, either through ordinary branching (lines 9-11) or the new one (lines 13-14). In either case, the number of nodes can increase only nitely-many times, since the branching causes the nodes to have (at most) one less integer rst-stage solution than their parent node. As a result, a xed rst-stage solution must be obtained in nite time, at which point, the node algorithm invokes the MICP solver instead of the PH algorithm. Due to A10, in nite time, the MICP solver either generates the next best solution for problemP-2, or proves that no better solution exists in the optimized node (perhaps, by showing CHAPTER 3. THE PH-BAB FRAMEWORK 41 Algorithm 5: A PH-BAB Algorithm for 2-Stage SMICPs with Pure Binary First-Stage Variables 1 Initialization: Create the root node o, letH fog. Set zU = +1, zL =1, ^ XU = NULL. 2 whileH6=; do 3 Pick node fromH and letH Hnfg. 4 Invoke Algorithm 4 to process . 5 if is feasible then // else, is fathomed 6 Get solution ^ X of , and a lower bound zL(). 7 if zL()<zU then // else, is fathomed 8 if ^ X has a fractional rst-stage solution then 9 Pick a fractional rst-stage component x j of ^ X . 10 Create nodes + and . Let Q + =Q =Q , while introducing the additional restrictions lstj =d^ x stj e;8s2FFt to Q + , and ustj =b^ x stj c;8s2FFt to Q . 11 LetH H[f ; + g. 12 else if ^ X has a fractional second-stage solution then 13 Create nodes f according to (3.11), and : according to (3.13). 14 LetH H[f f ; : g. 15 else 16 Obtain an upper bound zU () and the corresponding solution ^ X U . 17 if zU ()<zU then 18 Let ^ XU = ^ X U , zU =zU (). 19 Fathom:H Hnf 0 jzL( 0 )zUg 20 end 21 LetH H[fg. // is replaced in H 22 end 23 end 24 end 25 Let zL = minfzL()j2Hg. 26 end that it is infeasible). Since both of these are assumed to happen in nite numbers of steps, the nite convergence of Algorithm 5 follows. Remark 3.4. The results presented in this section can be extended to a special class of multi-stage SMICPs, where, at each stage, all nonanticipative decisions (that connect the stages) are pure binaries. The niteness of a corresponding algorithm can be shown with a similar argument to Proposition 3.5. 3.4. Computational Experiments In this section, we study the performance of the PH-BAB framework with comparisons to two other solution methods for SMIPs, which are commonly used in the literature. We will consider two PH-BAB implementations, where the rst one is an exact PH-BAB method, whereas the latter is an approximation. Our main studies will be conducted on instances which involve two stages and satisfy the pure binary rst-stage variables assumption of §3.3. Our choice of instances was motivated by the fact that, for this class of problems, benchmark instances are available and have been extensively studied in the literature. We begin with presenting our experimental design, which includes a discussion on the test instances, analyzed algorithms, and considered performance metrics. The rest of the section is devoted to the presentation of the outcomes of our computational experiments. 3.4.1. Experimental Design Our main experiments are conducted on two sets of stochastic server location problems (SSLPs). The rst set, due to Ntaimo and Sen (2005), has been used as a common benchmark in the stochastic programming literature. These instances are labeled according to SSLP.n.m.S where n and m represents certain problem parameters, whereas CHAPTER 3. THE PH-BAB FRAMEWORK 42 S stands for the number of scenarios. Each problem has n binary variables in its rst stage; nm binary and n continuous variables in the second stage. These instances can be obtained from Ahmed et al. (2015). As our second dataset, we introduce the quadratic extensions of these instances. In particular, we consider a quadratic objective function in the second stage of the problem. We label these instances as SSLP.Q.n.m.S. Further details on these instances can be found in Appendix B. We note that both of these datasets have pure binary variables in the rst stage, which allows the use of the developments given in §3.3. A second tier of experiments are also conducted on instances that do not satisfy the assumptions of §3.3, in order to demonstrate the potential of our framework on dierent families of SMICPs. To this end, we consider the ve network-design problem instances that appeared in Barnett et al. (2017), which are originally based on Ruszczy nski (2002). Each instance contains 30 binary and 30 continuous variables in its rst stage; 30 binary and 3,000 continuous variables in the second-stage. The number of scenarios in these instances are limited to 50, but the scenario subproblems are much more challenging due to their sizes, and the many disjunctions { modeled using Big-M's { that the formulations contain. We assess the performance of four dierent algorithms. The rst two of these are variants of PH-BAB imple- mentations which are specialized to two-stage problems with pure binary rst-stage variables (see §3.3). The third algorithm is the PH-Heuristic which is described in §3.1. Finally, we solve problemP-2 without any decomposition, through a deterministic mixed-integer programming solver. All algorithms, and the details within, are summarized in Table 3.1. Algorithm Description PH-BAB Algorithm 5. A single PH iteration is allowed per node. The default choice for is 100. PH-BAB-Apx An approximate variant of Algorithm 5. The PH algorithm optimizes nodal relaxations until a specied convergence tolerance is reached. No bound computations are performed on the nodal relaxations, and the nal objective value from the PH algorithm is assumed to be the optimal objective value of the processed node. The considered convergence criterion for the PH algorithm is P S s=1 pskxs xsk2 10 1 . The default choice for is 100. PH-Heuristic The PH algorithm which is executed with modied subproblems that are modeled as MICPs (see §3.1). The considered convergence criterion is P S s=1 pskxs xsk2 10 1 . Cycle detection is performed through monitoring the objective values from the last 10 iterations, and terminating the algorithm when all of them have already been recorded in earlier iterations. Unless specied otherwise, is selected adaptively, according to the sep heuristic given in Watson and Woodru (2010). Upon termination, an integral rst-stage solution is ensured by rounding the produced solution x to the nearest integer. This solution is used to compute an upper bound according to Proposition 3.3, with the integrality restrictions imposed. A lower bound is computed according to Proposition 3.2, again, with the integrality restrictions imposed (see Gade et al. (2016)). The computational time associated with these bound computations are not re ected into our reports. CPLEX Optimization of problemP-2, without decomposition, but with all of the default deterministic mixed-integer programming technologies of the solver CPLEX. Table 3.1: Descriptions of the tested algorithms. In all our experiments, the (relative) optimality gap is calculated according to (^ z UB ^ z LB )=j^ z UB j, where ^ z LB is the nal lower bound and ^ z UB is the nal upper bound produced by the corresponding algorithm. We consider a relative optimality tolerance of 10 4 and 10 2 for the linear and the quadratic test instances, respectively. All other tolerances are consistent with that of major commercial MIP solvers. In order not to underestimate the performance of the PH-Heuristic, we compute the relative optimality gap of the produced upper bounds (^ z UB ) according to (^ z UB z ? )=jz ? j, where z ? denotes the optimal objective value of the corresponding instance. When z ? is not available (for instance, for some quadratic instances), the objective value of the best upper bound (z ? UB ), obtained from any of the exact algorithms, is used instead. Finally, noting that PH-BAB-Apx is not an exact algorithm, we measure the error in its produced solution according to Error = max ( max(^ z LB z ? UB ; 0); max(z ? LB ^ z UB ; 0) ): CHAPTER 3. THE PH-BAB FRAMEWORK 43 Above,z ? LB denotes the best lower bound obtained from any of the exact algorithms. In simple terms, this measure indicates the maximum absolute deviation of the produced lower and upper bounds, from their (best available) theoretical upper and lower limits, respectively. Our runs were performed on a single thread of a Dell Desktop PC with Intel ® Core ™ i7-3770S CPU @ 3.10 GHz, 7.68 GB of RAM, and running Ubuntu Linux 12.04.3 LTS. For all algorithms, a time limit of 3600 seconds is imposed. The mathematical programs were solved through the Concert Technology class library of IBM ILOG CPLEX 12.6.2 (CPLEX). In the PH-BAB, PH-BAB-Apx, and PH-Heuristic algorithms, the linear and quadratic problems were solved using the dual simplex algorithm due to its superior performance. The remaining parameters of CPLEX were preserved at their default values. 3.4.2. Performance Comparisons on Linear SMIP Instances In Table 3.2, we present the solution times and the optimality gaps reported by all tested algorithms. Our main takeaways are presented in bullets, with ensuing paragraphs that provide further details. PH-BAB PH-BAB-Apx PH-Heuristic CPLEX Instance Optimal Objective Time Gap Time Gap Error Time Gap UB-Gap Time Gap SSLP.5.25.50 -121.600 0.6 - 0.5 - - 5.1 6.7 - 1.3 - SSLP.5.25.100 -127.370 1.2 - 1.1 - - 11.8 6.6 - 3.9 - SSLP.5.50.100 -323.700 1.7 - 1.5 - - 3.9 1.7 - 6.9 - SSLP.5.50.500 -320.800 8.1 - 7.0 - - 14.7 1.5 - 181.3 - SSLP.5.50.1000 -320.148 16.0 - 13.8 - - 29.5 1.4 - 523.7 - SSLP.5.50.2000 -320.602 32.0 - 27.8 - - 59.0 1.4 - 3600.1 42.5 SSLP.10.50.50 -364.640 9.4 - 8.7 - - 40.8 3.9 - 198.8 - SSLP.10.50.100 -354.190 19.5 - 18.6 - - 89.0 4.9 - 1685.6 - SSLP.10.50.500 -349.136 86.1 - 80.8 - - 331.6 4.7 - 3601.7 0.2 SSLP.10.50.1000 -351.711 172.4 - 163.6 - - 614.8 4.4 - 3600.2 64.3 SSLP.10.50.2000 -347.262 333.1 - 309.6 - - 1192.4 4.5 - 3600.2 112.6 SSLP.15.45.5 -262.400 2.0 - 1.4 - - 5.2 3.4 0.5 3.4 - SSLP.15.45.10 -260.500 10.9 - 2.4 - - 33.7 5.8 - 1.0 - SSLP.15.45.15 -253.600 5.4 - 3.2 - - 3601.5 8.7 0.2 11.9 - Table 3.2: Performance of algorithms for a class of stochastic linear programs (relative optimality tolerance = 0.01%; solution times are given in seconds, gaps are given as percentages). : Time limit of 3600 seconds is exceeded. Among all the contenders, the PH-BAB and the PH-BAB-Apx algorithms perform signicantly faster for almost all of the linear instances. We observe that the combined solution time of these algorithms on the whole dataset is less than 25 minutes, whereas the contender algorithms may spend more time than this on particular instances. The recorded speed-ups with respect to the PH-Heuristic average to a factor of four, even after excluding the outlier instance SSLP.15.45.15. Aligned with our expectations, the performance of CPLEX is the poorest when the problems involve more than a handful of scenarios. In particular, we report four cases where the algorithm cannot certify optimality within the time allotted. The one exception is SSLP.10.45.10 where CPLEX performs the fastest among all of the contenders. However, this instance consists of very few scenarios, and for such instances, we should not have any claim for superiority for any of the decomposition-based algorithms. These outcomes should not come at a surprise when the steps of each algorithm are carefully scrutinized. For instance, the iterations of both the PH-BAB and PH-BAB-Apx algorithms mainly comprise solving series of quadratic optimization subproblems. In contrast, the PH-Heuristic deals with non-convex mixed-integer quadratic subproblems. As a consequence, despite the benet of decomposition, its performance is weakened by its expensive iterations. In the case of CPLEX, the poor performance can only be explained by noting that it does not take advantage of the stochastic-programming structure of the tested problems. It is interesting to note that, even though CPLEX operates on a linear program, it can still be beaten by methods that perform quadratic optimizations. CHAPTER 3. THE PH-BAB FRAMEWORK 44 For all linear instances, the PH-BAB and the PH-BAB-Apx algorithms have produced optimal solutions. The PH-Heuristic algorithm was also successful in identifying optimal solutions with minor exceptions. Across the PH-based methods, solution quality does not appear to be a signicant determinant of the performances of the algorithms. Indeed, the self-reported optimality gaps of the PH-Heuristic underestimates its performance. In the UB-GAP column of Table 3.2, we see that the gaps between the objective values of the optimal and the produced solutions are all zero, with only two minor exceptions. Finally, with CPLEX, instances with large numbers of scenarios may result into gaps exceeding 100%. With respect to the PH-BAB algorithm, the PH-BAB-Apx spends on average 30% less computational eort on the linear instances. In the PH-BAB algorithm, the bound computations on the nodes { which provide the exactness guarantee to the algorithm { leads to more computational eort than the gain achieved from the partial processing of the nodes. We observe that the PH-BAB-Apx algorithm is faster than the PH-BAB, without any exception. This indicates for the tested instances that small numbers of iterations are sucient to attain the specied convergence tolerance in the PH-BAB-Apx algorithm. Tighter tolerances will surely invalidate the conclusion drawn in here. For applications where the optimality of the produced solution is not the utmost concern, the PH-BAB-Apx can be used as an eective heuristic algorithm. While the risk of sub-optimality may not be eliminated, our experiments reveal that the PH-BAB-Apx algorithm produces optimal solutions fast and without any error. Moreover, in comparison to the PH-Heuristic, the PH-BAB- Apx is not prone to oscillating iterates, or, as is the case for SSLP.15.45.15, unexpectedly-slow convergence rates. The algorithm has minimal tuning requirements (mainly, and the convergence tolerance for PH optimizations). At least for the class of instances we have experimented with, PH-BAB-Apx algorithm appears to be an eective solution approach. We now focus on the behavior of the PH-BAB algorithm on an individual instance. Figure 3.3, illustrates the progress of the algorithm for our largest instance SSLP.10.50.2000. The gray-shaded areas visualize the absolute optimality gaps computed after (partial) optimizations of nodal relaxations with the PH algorithm. We observe that in general the PH-BAB algorithm spends a uniform amount of eort on the individual nodes with one important exception. After the rst 300 seconds, no optimality gap is reported for almost 30 seconds. The underlying reason for this is the mixed-integer programming optimizations of the scenario subproblems (see §3.3). This procedure { however expensive it may be { results into the rst (and the optimal) integer-feasible solution for the considered problem. Immediately after this, the PH-BAB algorithm attempts to perform mixed-integer programming optimizations for another node. In this case, however, the rst phase of the proposed node algorithm detects in short amount of time that a lower bound for this node exceeds the objective value of the best integer- feasible solution (see Algorithm 4 in §3.3). Figure 3.3 also reveals that the optimality gaps of the nodes tend to be smaller as the PH-BAB algorithm proceeds. More importantly, even if the gaps are not small (for instance, as in the last processed node), the nite iteration limit on the PH algorithm and the computed lower bounds prevent the algorithm to spend unnecessary eort on nodes with no promise. 3.4.3. Sensitivity of PH-based Algorithms to In the literature, it has been shown that the convergence behavior of the PH algorithm diers according to the magnitude of . For a class of convex problems, Mulvey and Vladimirou (1991) observe that high values of lead toX moving towards nonanticipativity fast, but slowing down later as X approaches the optimal solution ^ X . Mulvey and Vladimirou (1991) also reports that smaller values of typically imply a more gradual convergence. For the mixed-integer case, Watson and Woodru (2010) report large numbers of iterations for small values of , supporting the latter claim of Mulvey and Vladimirou (1991). In this section, we aim to understand the eects of on all PH-based methods under consideration. We conduct tests on three dierent settings: = 10; 100; and 1000. All runs were performed on linear instances. For each algorithm, performance measures of interest are reported in Table 3.3. CHAPTER 3. THE PH-BAB FRAMEWORK 45 Figure 3.3: Progress of the PH-BAB algorithm and the optimality gaps of the processed nodal relaxations (SSLP.10.50.2000; = 100; single PH iteration is allowed per node) The selection of has the least impact on the PH-BAB algorithm's performance, among the performances of all other assessed PH-based methods. By design, the partial exploration of the nodes and the embedded bounding mechanisms allow the PH-BAB algo- rithm not to waste eort on nodes with little promise. For both small and large values of , the PH-BAB algorithm improves its precision, only when necessary, through increasing the number of nodes explored. In comparison, the eect of can be signicant for the PH-BAB-Apx and PH-Heuristic algorithms. For the PH-BAB-Apx algorithm, we observe that larger values of can signicantly speed up the convergence, but have the potential to reduce the accuracy. Table 3.3 reveals that the setting = 1000 leads to three cases with positive Error that amounts to nearly 2 units (< 1%, relative to the optimal objective values of the instances). This outcome can be explained with the ndings of Mulvey and Vladimirou (1991). When is large and the convergence tolerance is not suciently tight, rapid convergence of the iterates towards the non-anticipativity can lead to premature termination of the PH-BAB-Apx algorithm, although the iterates may still be far from the optimal solution. While the accuracy can be improved with a tighter convergence tolerance, this would necessitate more computational eort, and the tradeo between speed and accuracy could not be eliminated. Similar conclusions can be made for the PH-Heuristic as well. We observe that the solution times and the number of iterations are signicantly higher when is small, but the solutions produced by the method are always optimal. When is large, the solution times become very competitive with respect to the PH-BAB algorithm, however the optimality gaps of the produced solutions are far from satisfactory. 3.4.4. Performance Comparisons on Quadratic SMIP Instances In our experiments on quadratic instances, we abandon the adaptive selection scheme sep of the PH-Heuristic algorithm, as it was neither prescribed for a quadratic setting nor did it provide a better performance. We consider a weaker relative optimality tolerance, namely 1%, as tighter tolerances render a portion of the mixed-integer programs unsolvable with reasonable computational eort. The optimal objective values of the test instances are obtained using the PH-BAB algorithm, without imposing a time limit, and with further instance-specic tunings. For three instances, however, an optimal solution cannot be certied as the memory requirements exceeded the limit in our computing environment. We summarize the performances of all algorithms in Table 3.4. Even though the specied optimality tolerance is weaker, in general, we observe that the optimization of the quadratic instances require much more computational eort than their linear counterparts. The algorithms behave similarly as they did with the linear instances, however, their dierences are more pronounced in the quadratic setting. CHAPTER 3. THE PH-BAB FRAMEWORK 46 PH-BAB PH-BAB-Apx PH-Heuristic Instance Nodes PH-Itr Time Gap Nodes PH-Itr Time Gap Error Itr Time Gap UB-Gap = 10 SSLP.5.25.50 31 26 0.5 - 35 130 1.4 - - 14 7.8 2.5 - SSLP.5.25.100 27 24 0.8 - 35 96 2.2 - - 12 15.0 2.3 - SSLP.5.50.100 35 30 1.4 - 35 76 2.6 - - 4 5.0 0.5 - SSLP.5.50.500 35 30 6.9 - 35 82 13.5 - - 3 19.4 0.5 - SSLP.5.50.1000 35 30 13.1 - 35 81 26.5 - - 3 38.4 0.4 - SSLP.5.50.2000 35 30 27.2 - 35 82 53.5 - - 3 75.5 0.4 - SSLP.10.50.50 275 264 11.6 - 229 1011 31.3 - - 13 91.3 2.3 - SSLP.10.50.100 287 274 25.1 - 235 1122 70.2 - - 19 284.6 2.3 - SSLP.10.50.500 273 262 110.1 - 211 1035 310.3 - - 11 751.9 2.4 - SSLP.10.50.1000 275 266 211.2 - 217 1053 610.4 - - 8 1119.7 2.3 - SSLP.10.50.2000 255 250 370.2 - 199 1003 1170.2 - - 8 2147.5 2.2 - SSLP.15.45.5 263 252 2.7 - 201 1115 6.5 - - 7 4.0 1.1 - SSLP.15.45.10 467 455 14.1 - 171 1021 13.1 - - 76 256.8 3.4 - SSLP.15.45.15 495 486 17.1 - 113 592 13.7 - - 111 785.6 3.3 - = 100 SSLP.5.25.50 33 28 0.6 - 31 36 0.5 - - 24 4.3 12.9 2.2 SSLP.5.25.100 34 28 1.2 - 35 35 1.1 - - 14 6.7 10.1 1.4 SSLP.5.50.100 37 30 1.7 - 35 30 1.5 - - 5 3.2 2.7 - SSLP.5.50.500 37 30 8.1 - 35 31 7.0 - - 4 13.2 2.1 - SSLP.5.50.1000 37 30 16.0 - 35 31 13.8 - - 4 27.1 2.0 - SSLP.5.50.2000 37 30 32.0 - 35 31 27.8 - - 4 54.4 2.0 - SSLP.10.50.50 221 218 9.4 - 199 266 8.7 - - 12 30.4 3.9 - SSLP.10.50.100 225 220 19.5 - 213 286 18.6 - - 11 51.5 4.9 - SSLP.10.50.500 209 208 86.1 - 195 269 80.8 - - 7 188.5 4.7 - SSLP.10.50.1000 215 214 172.4 - 201 268 163.6 - - 8 412.2 4.4 - SSLP.10.50.2000 201 200 333.1 - 187 261 309.6 - - 10 918.9 4.5 - SSLP.15.45.5 181 176 2.0 - 153 186 1.4 - - 3 1.2 4.4 1.2 SSLP.15.45.10 205 203 10.9 - 141 213 2.4 - - 7 5.3 5.8 - SSLP.15.45.15 137 135 5.4 - 97 143 3.2 - - 7 11.6 8.6 - = 1000 SSLP.5.25.50 35 28 0.7 - 29 23 0.5 - - 4 0.8 12.9 2.2 SSLP.5.25.100 37 28 1.5 - 30 26 0.8 - 1.8 3 1.3 23.0 11.7 SSLP.5.50.100 37 30 1.9 - 33 29 1.3 - - 4 2.1 15.0 10.6 SSLP.5.50.500 35 30 8.1 - 33 29 6.1 - - 3 8.8 6.7 4.3 SSLP.5.50.1000 37 30 17.8 - 33 29 12.0 - - 3 17.9 6.8 4.5 SSLP.5.50.2000 37 30 35.8 - 33 29 23.9 - - 3 36.0 6.9 4.6 SSLP.10.50.50 259 252 11.7 - 143 141 4.7 - - 4 10.1 9.8 5.4 SSLP.10.50.100 265 258 22.6 - 153 149 10.0 - - 3 14.7 18.9 11.8 SSLP.10.50.500 239 234 102.6 - 125 125 36.8 - - 5 86.8 11.2 5.8 SSLP.10.50.1000 243 238 202.3 - 149 149 83.7 - - 4 176.9 10.7 5.7 SSLP.10.50.2000 229 228 370.8 - 125 125 145.6 - - 4 329.5 11.3 6.1 SSLP.15.45.5 181 179 2.9 - 43 45 0.4 - 1.2 3 0.9 4.4 1.2 SSLP.15.45.10 275 273 4.5 - 105 104 2.2 - 1.2 3 3.7 7.5 1.6 SSLP.15.45.15 185 181 10.3 - 61 60 1.9 - - 7 9.0 9.1 0.5 Table 3.3: Performance of the algorithms under three dierent settings (relative optimality tolerance = 1%; solution times are given in seconds, gaps are given as percentages). Among all the contenders, the PH-BAB and the PH-BAB-Apx algorithms perform signicantly faster for almost all of the quadratic instances. Both the PH-Heuristic and CPLEX may occasionally hit the alloted time limit and may produce low-quality solutions. A particular example is the SSLP.Q.10.50.2000 instance, where the PH-Heuristic hits the time limit, and produces a solution that has 14.8% worse objective value than the best-known upper bound. For this instance, the PH-BAB and the PH-BAB-Apx algorithms can provide a solution (within the desired tolerances) in half an hour. CPLEX, on the other hand, spends more time than this and terminates with no feasible solution. Similar errors have been observed for three other instances, which emphasizes the importance of decomposition on stochastic programming problems. CHAPTER 3. THE PH-BAB FRAMEWORK 47 PH-BAB PH-BAB-Apx PH-Heuristic CPLEX Instance Optimal Objective Time Gap Time Gap Error Time Gap UB-Gap Time Gap SSLP.Q.5.25.50 -117.160 0.7 0.7 0.6 0.7 - 3.4 6.6 0.0 3.2 0.9 SSLP.Q.5.25.100 -123.395 1.9 0.7 1.2 0.7 - 6.2 5.3 0.0 8.6 0.9 SSLP.Q.5.50.100 -319.025 3.6 0.8 2.7 0.8 - 13.3 3.1 0.1 178.3 1.0 SSLP.Q.5.50.500 -315.139 19.8 0.9 13.3 0.9 - 60.7 3.2 0.1 3600.1 4.2 SSLP.Q.5.50.1000 -314.936 39.9 0.8 31.2 0.8 - 122.1 2.9 0.1 1074.3 error SSLP.Q.5.50.2000 -315.171 79.5 0.7 62.5 0.7 - 248.7 3.0 0.1 2909.3 error SSLP.Q.10.50.50 -319.010 41.4 1.0 19.5 1.0 - 253.0 8.2 1.4 3600.0 1.4 SSLP.Q.10.50.100 -310.840 40.0 1.0 58.4 1.0 - 389.6 7.3 0.2 3600.1 4.0 SSLP.Q.10.50.500 -307.177 y 322.1 1.0 280.8 0.6 - 1702.3 6.8 0.0 1196.1 error SSLP.Q.10.50.1000 -309.126 y 639.6 0.9 286.4 0.9 - 3684.3 6.7 0.0 3600.4 123.8 SSLP.Q.10.50.2000 -305.045 y 1275.3 1.0 1100.3 1.0 - 4164.0 25.3 14.8 2173.5 error SSLP.Q.15.45.5 -186.700 3.2 0.9 2.4 0.8 - 40.7 0.7 0.0 173.4 1.0 SSLP.Q.15.45.10 -182.650 291.4 0.9 289.9 0.8 - 517.4 6.9 0.1 3600.0 1.5 SSLP.Q.15.45.15 -177.367 1162.7 0.7 838.8 0.7 - 732.7 8.4 0.2 3600.0 3.4 Table 3.4: Solution times and percentage optimality gaps for quadratic test instances (relative optimality tolerance = 1%; solution times are given in seconds, gaps are given as percentages) : Time limit of 3600 seconds is exceeded. y : A certicate of optimality is not available. 3.4.5. Parallel Computations in the PH-BAB Framework The PH-BAB framework is well-suited for parallelization, as it involves many subroutines (e.g., the PH algorithm, or possible MIP solves) which can be executed in parallel with minor need for synchronization. In this section, we demonstrate the performance of the PH-BAB and PH-BAB-Apx algorithms in a multi-threaded environment. We only consider parallel computations for independent subproblems, and leave parallel-search of the branch-and-bound tree as a future research direction. In order to process meaningfully-large numbers of subproblems in parallel, we change our computing environment. To this end, we use 32 threads of a Dell Desktop PC with Intel ® Xeon ® E5-2630 v3 CPU @ 2.40 GHz, 32 GB of RAM, running Ubuntu Linux 12.04.3 LTS. Note that this computer should still be considered as an ordinary desktop machine. Indeed, the individual cores of this computer demonstrate inferior performance than the cores used for the earlier experiments. The results of our experiment are presented in Table 3.5. PH-BAB PH-BAB-Apx PH-BAB PH-BAB-Apx Instance Serial Parallel Serial Parallel Instance Serial Parallel Serial Parallel SSLP.5.25.50 1.09 0.45 1.01 0.37 SSLP.Q.5.25.50 0.84 0.36 0.57 0.26 SSLP.5.25.100 2.02 0.85 1.91 0.51 SSLP.Q.5.25.100 2.10 0.55 1.47 0.38 SSLP.5.50.100 2.59 0.67 2.15 0.48 SSLP.Q.5.50.100 3.46 0.87 3.01 0.70 SSLP.5.50.500 11.86 1.36 9.85 1.06 SSLP.Q.5.50.500 19.49 2.05 15.28 1.46 SSLP.5.50.1000 23.62 2.22 19.46 1.92 SSLP.Q.5.50.1000 37.49 3.37 34.62 2.73 SSLP.5.50.2000 45.85 3.84 39.84 3.29 SSLP.Q.5.50.2000 74.90 5.88 62.58 5.08 SSLP.10.50.50 12.72 5.82 12.16 3.42 SSLP.Q.10.50.50 65.51 64.02 32.49 17.25 SSLP.10.50.100 25.61 5.68 27.07 5.97 SSLP.Q.10.50.100 58.33 18.76 107.16 39.25 SSLP.10.50.500 114.53 11.04 115.58 11.05 SSLP.Q.10.50.500 650.80 106.19 727.58 67.66 SSLP.10.50.1000 225.70 19.43 230.61 18.66 SSLP.Q.10.50.1000 1418.04 137.29 618.21 46.79 SSLP.10.50.2000 420.31 32.47 433.33 33.41 SSLP.Q.10.50.2000 2586.16 377.07 2272.73 198.54 SSLP.15.45.5 3.04 4.04 1.75 3.52 SSLP.Q.15.45.5 5.31 4.71 7.61 3.06 SSLP.15.45.10 8.51 9.40 2.91 2.62 SSLP.Q.15.45.10 3600.00 3603.13 3600.00 3602.63 SSLP.15.45.15 4.78 18.81 4.29 2.31 SSLP.Q.15.45.15 3600.00 3603.90 3600.00 3603.71 Table 3.5: Solution times of PH-BAB and PH-BAB-Apx for a class of stochastic linear and quadratic programs (optimality tolerances are kept consistent with the previous experiments; solution times are given in seconds). : Time limit of 3600 seconds is exceeded. Our experiment reveals that parallelization can indeed fulll the promised prospect of superior computational performance in the PH-BAB framework. As expected, the highest computational gains are observed for instances with the largest numbers of scenarios. In fact, for instances with small numbers of scenarios (e.g., SSLP.15.45.15), CHAPTER 3. THE PH-BAB FRAMEWORK 48 parallel runs may lead to inferior computational times. A typical explanation for this is the additional computational overhead for managing parallelization. Furthermore, the reader should also be reminded that an abundance of threads will lead to more optimizations that are not warm-started by the solutions of earlier subproblems. Indeed, as the number of parallel threads in a computing platform grows, the impact of warm-starts on subproblem- optimizations diminish. Such a drawback might be addressed with better coordination in parallelization. Even so, the overall improvements in the solution times for the majority of our test instances appear to be quite signicant. We are pleased to make the following observation: Our experiments reveal the decade-long progress made in (stochastic) mixed-integer programming: Problems which could not have been solved to optimality in two hours (Ntaimo and Sen, 2005) can now be optimized in half a minute. 3.4.6. Performance of PH-BAB Beyond Specially-Structured Problems Our nal analysis demonstrates the performance of our framework on a family of instances that does not satisfy the assumptions of§3.3. To this end, we consider the instances that were experimented with in Barnett et al. (2017). The presence of continuous variables in the rst stage of these problems makes them particularly challenging, as the convergence of the PH algorithm can no longer be aided by a branch-and-bound process. Accordingly, exact evaluations of the second-stage value functions may no longer be an option, since continuous variables cannot be xed unlike their binary counterparts. As a result, the PH-BAB algorithm may need to process signicant (but nite) numbers of nodes until discreteness is resolved, after which, the PH algorithm can run its course to produce feasible solutions. In view of large numbers of nodes explored during a branch-and-bound process, it is not an uncommon phe- nomenon for branch-and-bound schemes to employ heuristics to discover feasible solutions. Along these lines, we periodically invoke a feasible-solution heuristic which resembles Algorithm 4. In particular, after a partial node optimization, we round and x all rst-stage integer variables, and solve the second-stage subproblems as mixed- integer programs, subject to the variable bounds of the node. Then, the solutions to all integer variables are xed, and the PH algorithm is iterated on these subproblems to ensure the nonanticipativity of the (rst-stage) continuous variables. The scenario subproblems in these experimented instances are more dicult than that of the SSLP instances. To surmount the challenge, we use the same computing platform that we considered in §3.4.5, with parallel compu- tations enabled. As mentioned in§3.4.1, the experimented instances in this section contain many Big-M constraints, which can commonly induce poor linear-programming relaxations. In order to give emphasis to the quality of our lower bounds, we set = 1 so that the PH algorithm can produce better lower-bounding sequences. In contrast, in the PH iterations of our feasible-solution heuristic, we set = 100 so that the convergence of primal variables occurs faster. We invoke our heuristic after every 25 (feasible) node optimizations, and label a solution as nonanticipative according to the convergence criterion, P s2S p s kx s xk 2 10 1 . For brevity, the focus of this subsection will only be on the performance of the PH-BAB algorithm. Our main conclusion is stated below. Within the tested set of problems, the PH-BAB algorithm achieves rapid convergence to high-quality feasible solutions, and thereafter shows a gradual progress towards closing the optimality gap. In Figure 3.4, we observe that the most signicant improvements in optimality gaps occur within the rst 10 minutes of the PH-BAB algorithm. While the algorithm continues to close the gaps afterwards, the paces of improvements slow down, perhaps because most rst-stage variables have already been branched on, and branching on second-stage variables makes less dierence. In Table 3.6, we present snapshots of the best upper and lower bounds that the PH-BAB algorithm obtained within three dierent allocations of solution time. The table shows that the feasible solutions provided by the PH-BAB algorithm are high-in-quality, and improves only marginally with increasing computational budget. In contrast, the best lower bounds exhibit very gradual progress towards closing the optimality gaps. Such an outcome CHAPTER 3. THE PH-BAB FRAMEWORK 49 time (sec) 0 1200 2400 3600 4800 6000 7200 Relative Optimality Gap 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 1ef50 2ef50 3ef50 4ef50 5ef50 Figure 3.4: Progress of the relative optimality gaps reported by the PH-BAB algorithm. is not particularly surprising, as nding high-quality feasible solutions, quickly, and expending signicant eort on improving them, is a typical behavior of many (commercial) branch-and-bound algorithms (see, for instance, Barnett et al. (2017), for a discussion on the behavior of CPLEX for these problems). For benchmarking purposes, we also consider the performance of the framework suggested in Barnett et al. (2017). The metrics, presented in Table 3.6, are taken directly from their paper, where the computations therein were performed on a signicantly superior computing environment 2 . We observe that the upper bounds obtained by the PH-BAB framework (within 600 seconds) are superior than those obtained by this contender. Similar outcomes can also be drawn by comparing the performance of the PH-Heuristic, presented in Barnett et al. (2017), to the results obtained here. CPLEX 3 Barnett et al. (2017) PH-BAB { Best UBs PH-BAB { Best LBs Instance Best UB Best UB (Time) 600s 1800s 7200s 600s 1800s 7200s 1ef50 158,653 162,163 (18,161s) 161,682 161,682 160,730 114,990 119,816 122,616 2ef50 151,060 156,211 (15,330s) 152,100 152,100 151,827 111,410 115,733 118,287 3ef50 161,466 165,733 (33,390s) 164,426 164,426 164,426 115,932 121,104 124,214 4ef50 153,854 157,229 (18,150s) 155,944 154,794 154,794 113,576 117,763 120,260 5ef50 150,401 152,686 (8,556s) 151,581 151,581 151,006 112,231 116,487 119,055 Table 3.6: Performance of PH-type algorithms that consider a branch-and-bound framework. 3.5. Notes In this chapter, we have developed a new decomposition framework for computing optimal solutions to multi- stage SMICPs. Our framework has been specialized and implemented for a challenging class of two-stage SMICPs. Computational experiments reveal the eectiveness of the proposed framework on a common mixed-integer stochas- tic linear programming benchmark, and its newly-introduced quadratic counterpart. Our experiments also reveal the signicant progress made in (stochastic) mixed-integer programming, over the last decade. The PH-BAB framework appears to be an important candidate for solving stochastic UC problems. The current literature, addressing these problems, often consider a handful of scenarios for representing random variables. Such small sets of scenarios typically cannot accommodate uncertainty, and may, in fact, produce inferior decisions. 2 The authors use a 96-core Intel Xeon workstation with 2.1 GHz processors and 1 TB of RAM, where 50 instances of the PH algorithm are executed in parallel, and each instance is allowed to use up to 2 cores. 3 The best objective values for these instances are taken from Barnett et al. (2017), where the authors use CPLEX to solveP, without decomposition, for 10,000 seconds within their computing environment. CHAPTER 3. THE PH-BAB FRAMEWORK 50 On the other hand, the PH-BAB framework can eciently handle large numbers of scenarios, and with further improvements on its branch-and-bound process, it may scale well to realistic stochastic UC instances. This, in turn, can lead to statistically-optimal solutions, and get us one step closer to mitigating renewable-uncertainty in power-system operations. Chapter 4 Taming the Duck: Can Stochastic Optimization Help? For decades, human-made climate change has been an alarming research subject for scientists, and a contentious topic of debate for politicians. On one side of the argument, people have pushed renewable-energy as a solution to climate change, whereas others have denied the very existence of global warming, and in some cases, impeded the progress of renewable-energy to preserve fossil-fuel industries. Yet, and perhaps not to anyone's surprise, scientists have been proven right for their projections of increasing global temperatures, reducing Arctic sea ice levels, and the ubiquity of natural-disasters, among others. The need for reducing greenhouse gas (GHG) emissions has become so apparent that the world's nations signed the Paris Agreement in 2015, and strive to keep their promises, despite the agreement's non-binding nature. Furthermore, the biggest obstacle against the growth of renewable resources, namely their economies of scale, has been lowered due to new technological advances, leading the way for marked additions of renewable-capacity in recent years (The Guardian, 2017). In spite of such good news, an engineer's perspective towards ever increasing share of renewables in power grids is not completely enthusiastic. The availability of renewable resources depends on nature, and harvesting solar and wind is not trivial, at least, in comparison to the extraction, storage, and burning of fossil-fuels. Accordingly, with increasing renewable integration, new challenges emerge in power systems that need to be addressed with caution, so that the eciency and reliability of operations can be maintained. The goal of this chapter is to examine these challenges, and oer specic suggestions, primarily, in the context of power-system operations planning, and more specically, planning under renewable uncertainty. The organization of this chapter is as follows. In §4.1, we elaborate the adversities of renewable-energy resources which cause concern for some engineering communities. This section will also elucidate our choice of title for this chapter. Next, in §4.2, we provide a modeling framework for experimenting with daily power-system operations in typical grids. Then, we use this framework to conduct experiments in §4.3 on a test grid. These experiments will reveal the challenges faced by the power industry, and assess the impact of algorithmic novelties to overcome the uncertain nature of renewable resources. Concluding remarks are given in §4.4. 4.1. The Quandary of Renewable Energy Integration Lawmakers throughout the U.S. have mandated that a signicant percentage of electricity supply should be derived from renewable resources | each state has set its own goal, with California being the most aggressive, The research in this chapter was conducted in collaboration with Harsha Gangammanavar, PhD, from the Department of Engineering Management, Information, and Systems, Southern Methodist University, Dallas, TX (e-mail: harsha@smu.edu). 51 CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 52 requiring 50% by 2030 (SB-350, 2015). State and local authorities (e.g., ISOs) have commissioned studies to assess operational considerations such as system reliability, market design, storage technologies and other devices. A recent study, commissioned by California ISO (CAISO), suggests that for renewable-penetration levels beyond 33%, one can expect a great deal of over-generation, while also facing the possibility of curtailment of renewable energy during daytime, or load-shedding during sunset (Olson et al., 2015). The reason underlying such challenges is two-fold. First, wind and solar generators are non-dispatchable; that is, their output cannot be adjusted based on market conditions without external storage facilities. Consequently, the time periods in a day where there is signicant solar and wind supply may not coincide with periods of sucient demand. Moreover, energy storage can still not be performed at a sucient scale to mitigate these imbalances. Second, dispatchable generators typically have limited ramping capabilities, and therefore may not be able to handle sudden drops in renewable output. Accordingly, a sucient number of generators must be kept operational to ensure that supply meets demand, seamlessly. A popular illustration of this issue is given by the \duck-chart" of CAISO (see Figure 4.1). This gure depicts the net-load (total load minus solar and wind generation) over a timespan of increased renewable integration. A surplus of solar energy during daytime leads to a dip in net-load, followed by a signicant upward ramp of production achieved through conventional generators. In a grid with limited storage capabilities, excess supply poses signicant challenges. Over-generation during daytime could lead to negative prices in the market, resulting in large shipments of energy to neighboring states (e.g., Arizona), while paying these states to accept the surplus at home (e.g., California; see Penn, 2017). On the other hand, substantial ramping requirements can push the loss- of-load probability to unacceptable levels, and may even cause load-shedding in certain areas, which jeopardizes system reliability and performance. Green grid reliability requires flexible resource capabilities To reliably operate in these conditions, the ISO requires flexible resources defined by their operating capabilities. These characteristics include the ability to perform the following functions: • sustain upward or downward ramp; • respond for a defined period of time; • change ramp directions quickly; • store energy or modify use; • react quickly and meet expected operating levels; • start with short notice from a zero or low-electricity operating level; • start and stop multiple times per day; and • accurately forecast operating capability. Reliability requires balancing supply and demand The net load curves represent the variable portion that ISO must meet in real time. To maintain reliability the ISO must continuously match the demand for electricity with supply on a second-by-second basis. Historically, the ISO directed conventional, controllable power plant units to move up or down with the instantaneous or variable demand. With the growing penetration of renewables on the grid, there are higher levels of non-controllable, variable generation resources. Because of that, the ISO must direct controllable resources to match both variable demand and variable supply. The net load curves best illustrate this variability. The net load is calculated by taking the forecasted load and subtracting the forecasted electricity production from variable generation resources, wind and solar. These curves capture the forecast variability. The daily net load curves capture one aspect of forecasted variability. There will also be variability intra-hour and day-to-day that must be managed. The ISO created curves for every day of the year from 2012 to 2020 to illustrate how the net load following need varies with changing grid conditions. Ramping flexibility The ISO needs a resource mix that can react quickly to adjust electricity production to meet the sharp changes in electricity net demand. Figure 1 shows a net load curve for the January 11 study day for years 2012 through 2020. This curve shows the megawatt (MW) amounts the ISO must follow on the y axis over the different hours of the day shown on the x axis. Four distinct ramp periods emerge. 34,000 32,000 30,000 28,000 26,000 24,000 22,000 20,000 18,000 0 Megawatts 12am 3am 4 6am 7 9am 1012pm 13pm 6pm 9pm Net load - January 11 Hour 2012 (actual) 2013 2014 2020 2017 2015 2016 2019 2018 start stop start stop Figure 1 www.caiso.com | 250 Outcropping Way, Folsom, CA 95630 | 916.351.4400 CommPR/2016 © 2016 California ISO California Independent System Operator 2 Figure 4.1: CAISO's duck chart, predicting four emerging ramping patterns with increased renewable integration (California ISO, 2016). A number of solutions have been proposed to accommodate the growing ramping requirements in the grid. For instance, CAISO and California Public Utilities Commission propose \ exible resource adequacy criteria and must oer obligations" which requires California utilities to procure sucient upward ramping capabilities (California ISO, 2018a). These fast-ramping reserves use fossil fuels (mainly, natural gas) and therefore undermine the original drive towards renewable resources, more specically, the reduction in GHG emissions. Curtailment of renewable output has been an important tool in handling excess supply (Olson et al., 2015). In fact, the amount of curtailed solar and wind energy has been rapidly rising with increasing levels of renewable integration (California ISO, 2018b, see, also, Penn, 2017). This loss of opportunity produces negative economic and environmental consequences (California ISO, 2017b), and may discourage the progress towards greater renewable-integration. Adding further to these, the share of rooftop solar panels | whose outputs are not controlled by the ISOs | is also growing and poses an additional risk (Penn, 2017). Other suggestions to mitigate such issues include investing in energy storage technologies, improving regional coordination for consuming surplus generation, enhancing demand-response CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 53 initiatives for achieving higher elasticity in load, all of which bear signicant promise for the days to come. As these complications play out within the ISO communities, researchers in some environmental science programs have published a well-cited study where the authors suggest that renewable energy can be used to power 100% of the electricity demand by 2050 (Jacobson et al., 2015). This study has received signicant criticism from a large group of engineers (Clack et al., 2017; Bistline and Blanford, 2016), because Jacobson et al. (2015) did not accommodate many engineering-considerations such as transmission capacity, reliability, and the like. Meanwhile, on the policy side, California's push towards 100% renewable-energy by 2045 (SB-100, 2018) has been facing sti resistance, in part, due to reliability considerations, among other reasons (Nikolewski, 2017). In any event, it is important to recognize that just like any other infrastructure, reliability, ecient system-operation, and environmental considerations require us to strike a delicate balance between multiple facets, such as generator-mix, grid-characteristics, natural resources, and of course, market economics and public policy. While we do not intent to be mired in such debates, a certain issue in the discussed quandary of renewable integration draws our attention. Solar and wind energy are not only non-dispatchable, but also subject to signicant variability and uncertainty, as the availability of these resources depends on hard-to-predict atmospheric conditions. Actual duck-charts, depicted in Figure 4.2, demonstrate the variable nature of these resources at an aggregate level. Note that the impact of variability can get more severe at a granular level. More importantly, the duck-charts that appear in the literature are consistently depicted as \complete" gures, whereas the shape of the curve is never known in advance, and its predictions may be far from perfect (see Figure 4.3 for an illustration of uncertainty in the grid). As a result, a decision-maker, performing economic dispatches of electricity, must plan ahead and determine the adequate amount of over-generation, renewable curtailment, and an optimal generation mix, under information uncertainty. (a) January 30, 2018 (b) March 21, 2018 Figure 4.2: Aggregate load and net-load gures demonstrating the variability in solar & wind supply (Source: CAISO - Today's Outlook. Retrieved from http://www.caiso.com/TodaysOutlook/Pages/default.aspx; net-load is depicted by the lower line) (a) January 31, 2018, 1:50pm. (b) January 31, 2018, end of day. Figure 4.3: Aggregate load and net-load gures demonstrating the uncertainty in solar & wind supply (Source: CAISO - Today's Outlook. Retrieved from http://www.caiso.com/TodaysOutlook/Pages/default.aspx; net-load is depicted by the lower line) Several measures have been taken to mitigate the volatility of solar and wind resources. The FERC has mandated that electricity dispatch must take place at 10-15 minute intervals, instead of the usual hourly interval, which had been the norm for fossil-fuel generators (Order No. 764, FERC, 2012). Such shorter dispatching intervals can help the system track changes in renewable generation. Modern ISOs adopted sophisticated prediction models CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 54 (e.g., neural network and regression models) to put together accurate renewable-supply and load forecasts (see, for instance, California ISO, 2017a; New York ISO, 2017a). During real-time operations, such forecasts are frequently updated using new (weather) data, so that improved commitment and generation decision are built based on them. Despite the attempts to capture the volatility in solar and wind supply, this chapter argues that a signicant amount of information obtained from prediction models are not well-harnessed by the decision models used in industry. Typically, prediction models produce the best estimates for the realizations of the randomness in the form of a single expected time-series per stochastic process (e.g., the expected day-long, hourly, supply-estimates for an individual wind turbine). These time series are then fed into decision models that optimize the commitment and generation decisions. The nature of these decision models is deterministic, that is, they consider the given time-series as the only possible realization of the uncertainty, and produce decisions that are optimal with respect to this input. However, just like a weather forecast, the quality of solar and wind forecasts quickly deteriorate as the prediction-horizon grows. As a result, decisions that are presumed as optimal, may actually perform very poorly, when the uncertainty is revealed in a manner disregarded by the decision model. For capturing both uncertainty and variability, advanced decision models, such as stochastic programs, can be adopted. These models consider a large pool of scenarios that represent the joint realizations of the stochastic processes. Given the presence of a statistical prediction model, simulating such scenarios is not dicult, as it boils down to sampling from the error distribution that is readily available within the prediction model. By considering multiple scenarios, the decision model would be able to hedge against uncertainty and variability, and produce decisions that are closer to being statistically optimal. We refer the reader to Sen and Deng (2018) for a detailed treatment of statistical-optimality. For many power-system operations planning problems, stochastic optimization has already been proposed as an appealing option. The reader may review§1.1.1 for references that incorporate uncertainty into UC formulations. In the context of economic dispatch (ED) problems, Lee and Baldick (2013) consider a stochastic ED formulation and recommend its use to accommodate uncertain wind output, whereas Gangammanavar et al. (2015) demonstrate that advanced stochastic optimization models can produce statistically-optimal decisions and more economic outcomes. In what follows, we investigate whether these methodologies are capable of oering similar promising outcomes in daily power-system operations. 4.2. Models for Power-System Operations Planning 4.2.1. Modeling Daily Operations Power systems are typically managed in two-phases. The rst phase involves day-ahead planning, where the demand and renewable-supplies are estimated, generation-bids are collected, and a day-ahead unit commitment (DA-UC) problem is solved to clear the market. The DA-UC problem produces a low-cost commitment schedule for all generating units, typically, over a day-long planning-horizon with hourly resolution. Additional capacities, such as operating reserves, are also procured in the day-ahead market (Doherty and O'Malley, 2005; Morales et al., 2009). The second-phase of operations addresses the real-time requirements of the grid. Throughout the day, ED problems are periodically solved to determine the generation-schedules of committed units. These problems are typically solved for much shorter time-horizons but with ner resolution. Additional needs, such as precautionary reliability measures, are also assessed and addressed during real-time operations. In practice, ISO operations may signicantly vary across dierent authorities, resulting in models with varying degrees of complexity. For instance, Bonneville Power Administration, which primarily produces hydro-electricity, performs bulk-hourly generation-scheduling, and has sucient range of load-following, regulation, and ramping capabilities for handling within-hour imbalances (Makarov et al., 2010). Advanced ISOs, such as New York ISO, consider additional UC problems with ner resolution, in order to make adaptive (de)commitment decisions in real-time (New York ISO, 2017a). The CAISO, an authority that is under stress due to the state's ambitious renewable goals, makes use of a multitude of UC and ED problems, both for day-ahead planning and real-time CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 55 operations. Each of these problems serve a dierent purpose, may be solved periodically throughout the day, and may have resolutions as high as ve minutes (California ISO, 2017a). This level of sophistication, among other things, allows the use of up-to-date renewable-supply estimations in real-time, and helps maintain the reliability and cost-eectiveness of the grid. To this end, we adopt a three-layered hierarchical framework to approximate semi-realistic power-system opera- tions. At rst, we solve a DA-UC problem to determine hourly commitment-schedules for the majority of generating units. The DA-UC problem uses the best-available load and renewable-supply forecasts in day-ahead. Then, during the operation day, we consider a series of short-term unit commitment (ST-UC) problems to be solved at every hour. These problems produce commitment-decisions for the remaining (fast-start) generators, and generation-decisions for all generating units, over a 3-hour time horizon with 15-minute resolution. The ST-UC problems give the grid the necessary exibility to accommodate sudden drops in the actual solar and wind outputs. Finally, every 15 min- utes, we solve an ED problem to determine real-time dispatch amounts over a 1-hour time horizon with 15-minute resolution. The timeline of our framework is illustrated in Figure 4.4. These three problems usually cover the main components of power-system operations planning. Notice that our framework does not involve an energy market, or ancillary services such as ramping reserves. We assume that the latter can serve any unmet demand that is revealed during our planning process, albeit at a higher cost. Our Framework for Modeling ISO Operations ‣ We will mimic the current state of decision-making processes in typical ISOs with deterministic UC and ED optimizations. • DDD: Deterministic DA-UC, Deterministic ST-UC, Deterministic ED ‣ We will, then, introduce stochasticity into the optimization models, and analyze the impacts. • DDS, SDD, SDS, … ‣ All models will consider a transmission network. DA-UC (-08:00): ST-UC (00:00): -08:00 00:00 06:00 03:00 09:00 00:00 ED (00:00): Problem Time-Horizon Problem Resolution Figure 4.4: Timeline for optimization problems in our framework. In general, most ISOs operate the framework in Figure 4.4, by solving the UC and ED problems in a deterministic fashion. In our experiments, we will use this setup as our main benchmark. Then, we will investigate the potential of using stochastic optimization models, by gradually substituting these deterministic problems with their stochastic counterparts. These experiments aim to reveal the extent at which such algorithmic enhancements can mitigate the quandary of renewable-energy integration. Capturing all components of ISO operations in an academic experiment can be a formidable task. In this chapter, we are primarily interested in the \planning" aspect of power-system operations. In this respect, issues regarding the \prediction" of demand and solar/wind outputs will not enjoy a deep emphasis, despite their evident importance. In the remainder of this section, we discuss the methodology and procedures used to operate the framework in Figure 4.4. 4.2.2. Optimizing Commitment and Generation Schedules Every day, local authorities must go through tremendous optimization challenges to manage the grid in a reliable and cost-eective way. Surmounting these challenges is dicult due to the size of the grid, physics of transmission, physical limitations of generators, and information uncertainty, among other reasons. Modeling the planning problems as mathematical programs often leads to large-scale MIPs (e.g., for UC problems), LPs or nonlinear programs (e.g., for ED problems) which, by themselves, pose signicant challenges to the optimization community. For instance, early approaches to address UC problems involved the use of heuristics or relaxation-based approaches, which oered limited modeling support. In recent decades, however, ISOs have adopted state-of-the-art MIP solvers, which has provided signicant modeling exibility and the prospect of optimality for the problems at hand. In fact, a FERC study puts the amount of potential savings to $100 millions per annum (O'Neill et al., 2011). Today, modern ISOs heavily rely on such solvers for addressing both UC and ED problems. Accordingly, we will CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 56 use them as the baseline approach in our experiments, and as the backbone of algorithms that perform stochastic optimization. Mathematical Programming for Power-System Operations We begin our discussion with a description of the modeling components that will form the mathematical programming formulations of the UC and ED problems. A nomenclature, involving the denitions of considered sets, variables, and parameters, is given in Table 4.1. Sets B: buses. L: transmission lines. G: generators. Gr : solar and wind generators. Gc: conventional generators (Gc = GnGr ). G j : generators that are located in bus j2B. T : time periods. Parameters G max g : generation capacity of g2Gc. G min g : minimum generation requirement for g2Gc. G avail gt : wind/solar availability for g2Gr , in t2T . G max g : ramp up limit for g2Gc. G min g : ramp down limit for g2Gc. UTg : minimum uptime requirement of g2Gc. DTg : minimum downtime requirement of g2Gc. B ij : susceptance of arc (i;j)2L. D jt : load in bus j2B, in period t2T . R jt : reserve requirement in bus j 2 B and period t2T . F max ij : maximum permitted ow through arc (i;j)2L. c s g : start up cost of g2G. c f g : no-load cost of g2G (i.e., the intercept of the cost curve). c p g : variable generation cost of g2G (i.e., the slope of the cost curve). max j : upper bound on the voltage-angle at bus j2B. min j : lower bound on the voltage-angle at bus j2B. o g : penalty for over-generation by g2Gc. c g : penalty for renewable curtailment in g2Gr . u j : penalty for unmet demand in bus j2B. Decision Variables sgt: 1 if g2G is turned on in t2T , 0 otherwise. xgt: 1 if g2G is operational in t2T , 0 otherwise. zgt: 1 if g2G is turned o in t2T , 0 otherwise. G + gt : generation amount of g2G, in t2T , which is consumed by the grid. G gt : over-generation amount by g2Gc, in t2T . G gt : renewable curtailment in g2Gr , in t2T . F ij;t : electricity ow through (i;j)2L, in t2T . jt : voltage angle at j2B, in t2T . D shed jt : amount of unmet load at j2B, in t2T . Table 4.1: Nomenclature for the mathematical formulations UC, ED(x), SUC, and SED(x). Generator commitment decisions are often modeled using three sets of binary variables (x gt , s gt , z gt ) that indicate whether g is operational, turned on, and turned o in period t, respectively (Garver, 1962). The following constraints model feasible commitment schedules: x gt x gt1 =s gt z gt ; 8g2G; t2T; (4.1) t1 X j=tUTg +1 s gt x gt ; 8g2G; t2T; (4.2a) t X j=tDTg s gt 1x gt ; 8g2G; t2T: (4.2b) Constraint (4.1) forms the relation among the binary variables, whereas (4.2) (i.e., the turn on/o inequalities of Rajan and Takriti, 2005) model the minimum uptime/downtime requirements of the generators. Two sets of variables are used to model generation decisions. TheG + gt variable denotes the amount of electricity produced by generator g in period t, and consumed by the grid. The second variable G gt can assume two dierent meanings, depending on the type of generatorg. For solar and wind generators (g2G r ), these variables capture the amount of supply that is curtailed, whereas for all other generators (g2G c ), they represent the amount of electricity that is over-generated, in period t. We assume here that there exists a mechanism to consume over-generation at the buses where conventional generators are located. In more realistic settings, the over-generated electricity should be accounted for at certain locations where a consumer (e.g., a neighboring grid or energy-storage facilities) exists. Such information is not readily-available in the dataset we consider for our experiments. CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 57 All conventional generators (including hydro) must obey certain physical requirements for attaining feasible production-schedules. We take into account (i) generator capacities, (ii) minimum generation requirements, and (iii) ramping requirements: G min g x gt G gt G max g x gt ; 8g2G c ; t2T; (4.3) G min g G gt G gt1 G max g ; 8g2G c ; t2T: (4.4) Above (and in the ensuing discussion)G gt is used to simplify our exposition, whereG gt =G + gt +G gt ;8g2G;t2T . In the literature, the ramping constraints (4.4) are usually strengthened with binary variables to enhance the computational performance of MIP solvers. In our study, we incorporated some of the developments made in Damc-Kurt et al. (2015) and Chapter 2. For the purpose of conciseness, we do not present them in here. For solar and wind generators, the forecast and actual supply time-series are assumed to be capturing the physical requirements that these generators are subject to. Accordingly, we only need to impose a solar/wind availability constraint, which is given as follows: G gt =G avail gt x gt ; 8g2G r ; t2T: (4.5) The above constraint implies that the amount of available renewable capacity that is not consumed by the grid is considered as curtailment. Electricity transmission is modeled using three sets of variables which represent the electricity ow (F ij;t ), bus voltage angles ( jt ), and the amount of unmet demand (D shed jt ). We begin with the ow-balance equations: X i2B:(i;j)2L F ij;t X i2B:(j;i)2L F ji;t + X g2Gj G used gt +D shed jt =D jt +R jt ; j2B; t2T: (4.6) The above constraint also involves reserve considerations (R jt ), which are modeled as the sum of contingency and regulation requirements. We consider linear (direct-current) approximations of power- ows in our models. These are given, in terms of the bus voltage-angles, as follows: F ij;t =B ij ( it jt ); 8(i;j)2L;t2T; (4.7) min j jt max j ; 8j2B; t2T: (4.8) Finally, transmission capacities are imposed with the following constraints: F min ij F ij;t F max ij ; 8(i;j)2L; t2T: (4.9) We note that further modeling details that pertain to improving computational performance (such as symmetry breaking constraints or valid inequalities) are omitted in this presentation. Models and Algorithms for Deterministic Optimization Using the modeling components described above, we state our UC formulation as follows: min X t2T X g2G c s g s gt +c f g x gt +c v g G + gt + X j2B X g2Gj\Gc o g G gt + X g2Gj\Gr c g G gt + u j D shed jt ! (UC) subject to: (4.1) (4.9); (x gt ;s gt ;z gt )2f0; 1g 3 ; (G + gt ;G gt )2R 2 + ; 8g2G; t2T; F ij;t 2R; 8(i;j)2L; t2T; jt 2R; D shed jt 2R + ; 8j2B; t2T: CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 58 The above MIP formulation serves as the foundation of both DA-UC and ST-UC problem formulations, but mod- ications have been made to accommodate their dierences. For instance, the DA-UC problem does not commit certain fast-response generators, which can be ensured by setting x gt = 0 for such generators. Likewise, the ST-UC problem does not alter the commitment decisions made by the DA-UC problem, in which case, the x gt variables can be xed according to the solution of the DA-UC problem. Next, we present our ED formulation: min X t2T X g2G c v g G + gt + X j2B X g2Gj\Gc o g G gt + X g2Gj\Gr c g G gt + u j D shed jt ! (ED(x)) subject to: (4.3) (4.9); (G + gt ;G gt )2R 2 + ; 8g2G; t2T; F ij;t 2R; 8(i;j)2L; t2T; jt 2R; D shed jt 2R + ; 8j2B; t2T: Notice that the above LP problem is dened as a function ofx (i.e.,x gt ;8g2G; t2T , in vector form). Accordingly, the values of these variables must be xed in (4.3) and (4.5), prior to solving an ED the problem. Note that the constraints for idle generators (x gt = 0) need not be included into the ED formulation. We optimize the above formulations using a commercial mathematical programming solver (i.e., CPLEX), which can handle deterministic MIPs and LPs with remarkable computational performance. Models and Algorithms for Stochastic Optimization Our stochastic optimization models will retain the same modeling components as their deterministic counterparts, but will include solar/wind availability as a random input. Following the notation in Chapter 1, we consider the random vector , whose realization under scenario s is s . The vector s has s gt ;8g2G r ; t2T , as its components, where s gt denotes the realization of solar/wind availability under scenario s, for generator g, in period t. For both UC and ED problems, we consider a two-stage stochastic programming framework to perform stochastic modeling and optimization. The rst stage of the stochastic UC (SUC) problems involves generator-commitment decisions, whereas the second-stage problem makes economic-dispatch decisions. Based on this description, our SUC formulation is given as follows: min X t2T X g2G c s g s gt +c f g x gt +E ED(x;) (SUC) subject to: (4.1) (4:2); (x gt ;s gt ;z gt )2f0; 1g 3 ; 8g2G; t2T; where, ED(x; s ) = min X t2T X g2G c v g G + gt + X j2B X g2Gj\Gc o g G gt + X g2Gj\Gr c g G gt + u j D shed jt ! subject to: (4.3); (4.4); (4.6) (4.9); G gt = s gt x gt ; 8g2G r ; t2T; (4.10) (G + gt ;G gt )2R 2 + ; 8g2G; t2T; F ij;t 2R; 8(i;j)2L; t2T; jt 2R; D shed jt 2R + ; 8j2B; t2T: Note that the second-stage problem (ED(x; s )) has the same form as theED(x) problem, except for the additional input s and the corresponding solar/wind availability constraint (4.10). The input x in ED(x; s ) only aects constraints (4.3) and (4.10) of this problem. CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 59 For stochastic ED (SED) problems, all decisions associated with the initial time period are made in the rst stage, whereas the decisions corresponding to latter periods are made in the second stage. Accordingly, our SED formulation is stated as follows: min X g2G c v g G + g1 + X j2B X g2Gj\Gc o g G g1 + X g2Gj\Gr c g G g1 + u j D shed j1 +E ED 0 (x;) (SED(x)) subject to: (4.3) (4.9); (excluding t = 1); (G + g1 ;G g1 )2R 2 + ; 8g2G; F ij;1 2R; 8(i;j)2L; j1 2R; D shed j1 2R + ; 8j2B; where, ED 0 (x; s ) = min X t2Tnf1g X g2G c v g G + gt + X j2B X g2Gj\Gc o g G gt + X g2Gj\Gr c g G gt + u j D shed jt ! subject to: (4.3); (4.4); (4.6) (4.9); (excluding t = 1); G gt = s gt x gt ; 8g2G r ; t2Tnf1g; (4.11) (G + gt ;G gt )2R 2 + ; 8g2G; t2Tnf1g; F ij;t 2R; 8(i;j)2L; t2Tnf1g; jt 2R; D shed jt 2R + ; 8j2B; t2Tnf1g: Similar to ED(x; s ), the input x in ED 0 (x; s ) alters the constraints (4.3) and (4.11). Note that the components ofx corresponding tot = 1 is not used in theED 0 (x; s ) problem, but preserved as an input for ease of exposition. We consider two dierent optimization algorithms that are applicable to the respective natures of the above mathematical programs. The SUC formulation is a SMILP, where integer variables only appear in the rst stage of the problem. Given a xed number of scenarios, such problems can be eciently handled by the L-shaped algorithm of Van Slyke and Wets (1969) (see, also, Wollmer, 1980, or in the context of deterministic optimization, Benders' decomposition (Benders, 1962)). The SED problem, on the other hand, is an SLP, which is somewhat less notorious than SMILPs due to its convex nature. This motivates us to adopt sampling-based algorithms which are guaranteed to converge to statistically-optimal solutions. To this end, we optimize SEDs with the stochastic decomposition algorithm of Higle and Sen (1994) (see, also, Sen and Liu, 2016). 4.2.3. Predicting Solar and Wind Outputs The accuracy of solar and wind forecasts depends on the quality and abundance of information collected from operation sites, and the ability of the prediction models to harness them. Modern ISOs already receive reliable streams of data from solar and wind farms, and new mechanisms are being deployed to monitor behind-the-meter generation (e.g., by rooftop solar panels; see, for instance, New York ISO, 2017b). The prediction models may directly process physical input (e.g., weather, satellite images) and convert it into power forecasts, use real-time and historical supply information to make statistical predictions, or do a combination of both (Monteiro et al., 2009; Tuohy et al., 2015). The produced forecasts may have a lookahead ranging from a few minutes to several days, depending on the nature and purpose of the forecasting model. For the purposes of our study, we assume that day-ahead supply forecasts for each wind and solar generator are available as input. However, we still need a scheme to model the potential deviations from these predicted supply time-series. This necessity can be addressed with an m-th order vector autoregression (VAR) model, which can be CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 60 stated as VAR(m) = y t+i = m X j=1 j y tj+1 + t+i ; 8t; i = 1; 2;::: where, y t is the day-ahead supply forecast for a vector of generators, t is a white-noise process with a time- invariant positive-denite covariance matrix; and ( j ) m j=1 are autoregression coecient matrices. In contrast to a univariate autoregressive process, a VAR model can capture both temporal and spatial correlation of the stochastic processes which govern renewable output. We refer the reader to Lutkepohl (2007) for further details on multivariate autoregressive processes, Hering et al. (2015) for more advanced VAR models, and Gangammanavar and Sen (2017) for a similar application of VAR models to power-system planning problems. In our study, separate VAR models are used for wind and solar generators, both based on two weeks of normalized past data. For creating scenarios, we sample from the distribution of the white-noise process t , and add these sampled-deviations to the forecasts corresponding to the day ahead ( y t ). The procedure is repeated from scratch for every tested day. For alternative scenario generation procedures, the reader is referred to Rios et al. (2015); Woodru et al. (2018) for solar, and Morales et al. (2010) for wind outputs, among others. The accuracy of day-ahead supply forecasts and scenarios deteriorate as the prediction-horizon grows. To prevent unlikely forecasts and scenarios from becoming input to the optimization models, it is crucial to update these time-series based on the current state of generator outputs. In this study, we consider the following update rule: ~ y t+i = i y t + ( y t+i y t+i1 ) + (1 i ) y t+i ; 8t; i = 1; 2;::: where, ~ y t and y t are the updated supply forecast and actual supply for a vector of generators, respectively, and f i g is a sequence of coecients with i 2 [0; 1];8i, and i & 0. The above equation adds the forecasted increases ( y t+i y t+i1 ) to the most-recently observed generator outputs (y t ), and sets the updated forecast as a combination of these and the original day-ahead supply forecasts ( y t+i ). The latter operation ensures that the updated forecasts can reach to day-ahead forecast levels, although the predictions for the timing of production-ramps are inaccurate. Similar to solar and wind output, the demand at each bus can be consider as a stochastic process, and the above models can be used for their modeling. However, ISOs can often obtain very accurate load forecasts during the day, and the errors in day-ahead load forecasts are less severe than that of the renewable resources. Owing to this, we disregard the uncertainty associated with the loads in the ST-UC and ED problems. In the DA-UC problem, however, we consider the day-ahead forecasts of the loads. Finally, when necessary, we use spline-interpolation to obtain time-series with ner resolution. 4.2.4. Simulating Power-System Operations In the last part of the modeling discussion, we describe how the input data of §4.2.3 and problems in §4.2.2 are brought together to operate the experimental framework in §4.2.1. We begin with the case of deterministic DA-UC optimization. For each experimented day, we formulate the problem using the day-ahead supply- and load-forecasts. This problem determines the hourly commitment decisions for a (specied) majority of generators, where the decisions are assumed to be strictly enforced during real-time operations. A stochastic version of this problem will only require the substitution of the day-ahead supply-forecasts with the sampled supply-scenarios. During the day, we either use the actual renewable supplies or their updated forecasts, depending on the time-lag between the start of an optimization and the (real-time) execution of its decisions. To this end, the initial periods (t = 1) of all ST-UC and ED problems are formulated using the true realizations of solar and wind availability, whereas the remaining periods are formulated according to the updated forecasts. Given the quality of shorter-term load forecasts, we opted for the use of the actual demand values in all periods of these problems. Each ST-UC problem determines the subhourly commitment decisions of the fast-response generators. These CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 61 decisions are assumed to be unaltered for the next hour of operations, and thereafter, the subsequent ST-UC prob- lem's decisions are adopted. We do not impose any reserve requirements for the rst period of the ED problems. We assume that the rst-period generation decisions of the ED problems approximate the average real-time generation that will be dispatched in the next 15 minutes of operations. Accordingly, all reported statistics (e.g., costs) will be based on the (rst-period) generation-decisions of ED problems, along with the (xed) commitment decisions of the DA-UC and ST-UC problems. For the stochastic optimization of ST-UC and ED problems, we substitute all updated supply-forecasts with the updated supply-scenarios, following the similar case of the DA-UC problem. Note that the ow of the simulation is not altered by the choice of optimization paradigm for the individual problems. The initial operating statuses of generating units depend on the experimental setting that we consider. For the majority of our experiments, the DA-UC problems assume that all generators can generate arbitrary (but feasible) amounts of output, or be turned o, at the initial time period. This allows us to assess somewhat independently- formed operation plans in week-long experiments. We will also conduct experiments where we analyze the long-run eects of a particular variable in our experiment. In such cases, the commitment and generation decisions of the previous day are passed as an input to the subsequent DA-UC problem. Implementing a rolling-horizon framework, involving the interaction of multiple problem-types, comes with an important caveat. Varying problem resolutions and time-horizons, along with updated input parameters, may lead to divergence of decisions across the UC and ED problems, which may result into infeasibility during optimization (e.g., a generator may not be able to ramp down even though it is scheduled to turn o). In our framework, we detect such situations prior to each optimization, and change the input parameters to ensure a smooth transition across all solved problems. 4.3. Computational Experiments The goal of our experiments is two-fold. Using our experimental setup, we will demonstrate the impact of increasing renewable integration into power systems. In particular, we will explore how operating costs, GHG emis- sions, and reliability-related statistics (such as over-generation and unmet demand) are aected by varying degrees of solar and wind penetration. Second, we will investigate whether the observed outcomes prompt us to upgrade our decision-making mechanisms. More specically, we will compare deterministic and stochastic optimization paradigms for managing power-system operations, and elucidate their impacts using the aforementioned metrics. We conduct our experiments on a small power-system, referred to as NREL118, which was introduced by a group of engineers in the National Renewable Energy Lab (NREL) for experimenting with renewable-energy resources (Pena et al., 2017). The topology of this system is based on the IEEE118 dataset, which has been widely used by the academic community for conducting experiments on power networks. The NREL118 instance contains 327 generators (75 solar generators and 17 wind turbines), 118 buses, and 186 transmission lines, along with the forecasted and real-time outputs of solar and wind generators. We note that this power system should be considered as \progressive", as it involves rich solar, wind, and hydro resources, and lacks the traditional fossil- fueled generators, such as nuclear and coal (see Figure 4.5). For further details on the instance, we refer the reader to the original paper. We consider three major variables in the design of our experiments. These are (i) the amount of solar and wind resources within the grid, (ii) the extent of reserve considerations, and (iii) the choice of modeling and optimization for the UC or ED problems. We assess three levels of the rst variable, where the rst corresponds to the values in the original set (NREL118 S/W), whereas in the latter two, we multiply solar and wind outputs (and capacities) with 2.5 (Medium S/W) and 5 (High S/W). Figure 4.6 illustrates a breakdown of average supply-capacities under these scenarios. Note that this gure pertains to \capacities", and should not be confused with the requirements of renewable-integration goals, as these are concerned with \the percentage of demand fullled by renewable resources". We set the reserve amounts (R jt ) as percentages of the input demand (D jt ) at each problem. Similar to before, we consider three cases, where the reserve-percentages are set as 0% (N/A), 2.5% (Medium), and 5% (High). Finally, CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 62 Hydro 11% Natural gas 68% Oil 2% Other 1% Solar 14% Wind 4% Figure 4.5: Total generation capacities in NREL118 with respect to fuel-types. Hydro 10.1% Natural gas 82.3% Oil 1.8% Other 0.9% Solar 2.7% Wind 2.2% (a) NREL118 S/W Hydro 9.4% Natural gas 76.6% Oil 1.7% Other 0.9% Solar 6.3% Wind 5.1% (b) Medium S/W Hydro 8.5% Natural gas 68.8% Oil 1.5% Other 0.8% Solar 11.4% Wind 9.1% (c) High S/W Figure 4.6: Breakdown of average supply-capacities in NREL118 with respect to fuel types (based on February data). we adopt a three-letter convention to represent the modeling and optimization choices for DA-UC, ST-UC, and ED problems, where the letters (D or S) respectively indicate their (deterministic or stochastic) nature. For instance, an optimization setting with stochastic DA-UC, deterministic ST-UC, and stochastic ED problems, will be labeled as SDS. We begin our experiments with an illustration of the duck-charts based on the NREL118 dataset, under varying solar and wind integration levels, in Figure 4.7. We observe that all net-load curves manifest signicant upward- ramping requirements. We note that the Medium S/W scenario somewhat resembles the situation that California faces. On the other hand, the High S/W scenario is rather extreme, and, perhaps, represents the fretted situation of over-investment in solar and wind energy for reaching renewable goals (see, for a brief discussion on over-investment, California ISO, 2017b). Under this scenario, we even observe massive amounts of negative net-load, suggesting that renewable curtailment (or export) is unavoidable (see, for instance, day 7 in Figure 4.7). Given the physical laws governing electricity transmission, renewable curtailment has also been a necessity for the case of Medium S/W. In Table 4.2, we present the amounts of unmet demand that have emerged during our planning processes. The results suggest that increasing solar and wind integration may increase the risk of load shedding, unless the reserve considerations are also tightened in accordance with the growing shares of these renewables. Moreover, the required amounts of reserves to maintain tolerable unmet-demand levels appear to be growing very fast, and eventually, could lead to increasing operational costs and heavy reliance on ancillary generation services. For instance, under the DDD setting, we were able to reach tolerable unmet-demand levels, only, after the reserve-requirements were set as 10% of the demand, which is double the percentage of what we consider for the High-reserve case. The use of stochastic DA-UC problems could mitigate such concerns. We observe from Table 4.2 that the SDS setting signicantly reduces the unmet-demand levels to tolerable ranges. Note that in contrast to a xed reserve requirement, a stochastic DA-UC problem can dynamically adjust its decisions based on the amount of uncertainty that lies in the day-ahead. Accordingly, in conjunction with reasonable reserve considerations, the SDS setting has the most potential to achieve reliable and economic operations. With improved scenario-generation processes, it may even be possible to eliminate unmet demand completely from our planning process. CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 63 -6000 -4000 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 Megawatts Load Net Load (NREL118 S/W) Net Load (Medium S/W) Net Load (High S/W) Figure 4.7: Aggregated loads and net loads (under three solar/wind integration scenarios) in the NREL118 dataset plotted across 7 days (based on February 14-20 data). NREL118 S/W Medium S/W High S/W Reserves DDD DDS SDS DDD DDS SDS DDD DDS SDS N/A 701.6 688.2 592.1 1,706.4 1,661.9 628.2 2,991.1 2,854.9 509.6 Medium 0.0 19.4 27.3 396.6 367.8 121.9 980.0 839.3 285.4 High 0.0 0.0 0.0 0.0 59.5 53.0 426.7 482.6 28.8 Table 4.2: Average daily unmet demand in terms of megawatts per day (based on February 14-20 data and experiments). As a second reliability metric, we report the over-generation amounts by conventional generators in Table 4.3. Our results conrm the outcomes of earlier studies, suggesting, increasing solar and wind integration may lead to more over-generation in the power system. We observe that the share of over-generated energy grows substantially with more renewable integration. On the other hand, we did not observe a consistent trend with respect to the modeling and optimization paradigms that we assess. Under the Medium S/W scenario, we observe that the DDS and SDS settings markedly reduce the amount of over-generation. In contrast, under the High S/W scenario, the same settings lead to higher over-generation amounts. Such a discrepancy may be explained by the dynamics of the grid. In particular, in our experiments, we do not consider any modications to the grid, such as the expansion of transmission capacities, which must, ideally, go hand in hand to accommodate the tested levels of renewable integration. NREL118 S/W Medium S/W High S/W Reserves DDD DDS SDS DDD DDS SDS DDD DDS SDS N/A 4.5 0.4 0.7 1,526.8 171.7 244.3 3,537.9 4,094.4 3,426.4 Medium 0.0 7.9 0.3 1,448.5 194.4 232.5 2,731.7 4,607.7 3,664.5 High 0.0 9.6 2.0 1,686.5 119.0 256.8 3,356.8 5,418.2 4,231.2 Table 4.3: Average daily over-generation by conventional generators in terms of megawatts per day (including hydro; based on February 14-20 data and experiments). In Table 4.4, we present the renewable curtailment statistics. As expected, renewable curtailment and renewable integration demonstrate a positive correlation. Especially in the case of High S/W, the curtailed amounts grow to be substantial, in part, due to the large share of solar resources in our dataset. This, exposes an impact of unbalanced investment in renewable-energy resources. We also observe that stochastic optimization models have the tendency to curtail more solar and wind energy. Such an outcome is not particularly surprising, as hedging against uncertainty should, in general, require the use of more reliable resources. We turn our attention to operating costs, which involve commitment and generation expenses, as well as cost- CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 64 NREL118 S/W Medium S/W High S/W Reserves DDD DDS SDS DDD DDS SDS DDD DDS SDS N/A 0.0 0.0 0.0 10,028.8 11,812.1 12,112.9 77,096.8 81,363.7 79,578.7 Medium 0.0 0.0 0.0 9,749.0 11,978.0 11,670.1 76,531.3 80,289.1 81,671.1 High 0.0 0.0 0.0 9,768.8 12,282.6 11,866.5 75,803.1 81,207.9 79,747.3 Table 4.4: Average daily solar and wind curtailment in terms of megawatts per day (based on February 14-20 data and experiments). estimates for managing under- and over-generation. We compute no-load and generation costs based on the complete cost-curves of generators. The start-up costs are simply added according to generators' state transitions. We report two types of costs that dier in the way over-generation and unmet-demand expenses are incorporated. To obtain an optimistic cost-estimate, we assume that all unmet demand is covered by natural-gas burning generators. Each megawatt hour (MWh) of unmet demand is charged with the highest variable-cost across all natural-gas generators in our dataset ($57.30). We ignore potential no-load costs for keeping such generators operational. Likewise, over- generation is assumed to be removed at no cost, perhaps, due to the availability of energy-storage facilities. In the pessimistic scenario, we consider the possibility of actual load-curtailment in the grid, which could amount to $50,000 per MWh (Olson et al., 2015). Furthermore, for over-generated electricity, we consider an export price of $25 per MWh, which represents the highest price that California has reportedly paid for eliminating excess generation (Penn, 2017). The nal cost gures are summarized in Table 4.5. NREL118 S/W Medium S/W High S/W Reserves DDD DDS SDS DDD DDS SDS DDD DDS SDS Optimistic N/A 10,875 10,876 10,874 9,480 9,510 9,510 7,842 7,910 7,927 Medium 10,942 10,944 10,948 9,483 9,510 9,599 7,940 8,048 8,094 High 11,050 11,060 11,049 9,559 9,626 9,752 8,017 8,160 8,302 Pessimistic N/A 19,636 19,468 18,266 30,796 30,262 17,355 45,210 43,580 14,312 Medium 10,942 11,186 11,289 14,444 14,104 11,122 20,193 18,556 11,680 High 11,050 11,060 11,049 9,570 10,370 10,415 13,365 14,219 8,688 Table 4.5: Average daily costs, given in terms of thousand dollars (based on February 14-20 data and experiments). We begin with the discussion of the optimistic cost-estimates. It appears that the operating costs decline with growing solar and wind integration, underlining the economic benets of renewable integration. Furthermore, we observe that tighter reserve requirements and the use of stochastic optimization models slightly increase operating costs. This outcome should not be misinterpreted by the reader, as potentially-damaging outcomes are not included in these \optimistic" estimates. In fact, the largest pessimistic cost-estimates are observed for cases where the reserve requirements are ignored and the DDD setting is adopted. These results further indicate that the adoption of stochastic optimization models can potentially lead to signicant gains in reliability, which could translate into economic savings, in the long run. Aside from the economic benets, the main drive towards renewable-energy resources comes from their environ- mental impacts, or to be accurate, the lack thereof. We assess these eects based on three types of GHG emissions and summarize our results in Table 4.6. To estimate these gures, we use the complete emission-curve of generators. Furthermore, we assume that any unmet demand is covered by a natural-gas generator with the highest emission factors. We ignore potential emissions of GHGs for keeping generators operational. With increased renewable integration, we see a clear trend of reduced emissions, especially in terms of CO 2 and NO x amounts. The results indicate that the use of stochastic models lead to slightly larger amounts of pollutants, as these models, in general, lead to more renewable curtailment. To demonstrate the daily variability in power-system operations, we focus on a single scenario with Medium S/W integration and Medium reserve-considerations. Table 4.7 summarizes the unmet demand and operating cost statistics per each day of the experiment. These results are based on a dierent set of experiments than the ones CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 65 NREL118 S/W Medium S/W High S/W Reserves DDD DDS SDS DDD DDS SDS DDD DDS SDS CO 2 N/A 105,232 105,221 105,237 89,833 90,057 90,565 70,812 71,424 73,498 Medium 105,697 105,688 105,876 90,089 90,374 91,530 71,543 72,472 74,796 High 106,281 106,275 106,518 90,618 91,204 92,582 72,259 73,539 76,150 NO x N/A 73,142 73,116 73,044 62,554 62,677 62,811 49,747 50,141 51,262 Medium 73,347 73,282 73,438 62,578 62,713 63,408 50,237 50,822 52,029 High 73,560 73,468 73,590 62,788 63,102 63,975 50,433 51,268 52,845 SO 2 N/A 1,754 1,753 1,749 1,654 1,653 1,651 1,553 1,554 1,556 Medium 1,751 1,747 1,750 1,653 1,651 1,653 1,563 1,564 1,550 High 1,737 1,732 1,737 1,647 1,645 1,649 1,448 1,452 1,545 Table 4.6: Average daily GHG emissions, given in terms of metric tonnes (CO 2 ) and kilograms (NOx, SO 2 ) (based on February 14-20 data and experiments). considered above. In particular, here, we pass generator histories across the days, therefore the solutions to the DA-UC problem solved in Days 2 to 7 may depend on the solution of the earlier problems. Table 4.7 illustrates the value of stochastic optimization models, where we observe increased reliability and reduced operating costs. Day 1 Day 2 Day 3 Day 4 Day 5 Day 6 Day 7 Average Unmet Demand DDD 45.6 683.9 706.6 0.0 0.0 189.9 1,604.8 461.5 DDS 55.4 139.0 1,224.6 147.0 0.0 154.1 1,189.1 415.6 SDS 0.0 266.5 269.9 80.0 0.0 238.8 893.3 249.8 Operating Costs (optimistic) DDD 9,195 9,089 9,179 10,468 10,500 9,836 9,014 9,612 DDS 9,192 9,147 9,442 10,372 10,489 9,855 8,933 9,633 SDS 9,350 9,471 9,224 10,512 10,629 10,019 8,938 9,735 Operating Costs (pessimistic) DDD 9,896 17,718 18,048 10,599 10,509 12,218 29,153 15,448 DDS 10,015 10,984 24,787 12,330 10,493 11,806 23,890 14,901 SDS 9,481 12,885 12,647 11,635 10,632 13,057 20,201 12,934 Table 4.7: Daily unmet demand and operating costs in terms of megawatts per day and thousand dollars, respectively (based on February 16 data and experiments; generator histories are carried across the days) Finally, we move our attention to the generation-mix which fullls the demand. We begin with an illustration of the potential of stochastic optimization for accommodating uncertainty. In Figure 4.8, we demonstrate the portions of demand fullled with respect to each fuel type, along with over-generation and renewable-curtailment amounts. We observe that the DDD setting cannot accommodate a sudden increase in renewable output during the day, and leads to over-generation. Towards sunset, the same setting cannot accommodate the ramping requirements in the presence of uncertain energy resources, and this time, leads to unmet demand. While increasing reserves aid the reliability of operations, the best outcome is observed when medium reserve-requirements and the SDS setting are adopted in tandem. In general, the percentages of demand fullled by renewable-resources appear to be almost perfectly correlated with the amount of renewable resources in the system. The dierences in curtailment amounts across dierent optimization paradigms only slightly eect these percentages, and never exceed a single percentage point. We observe, on average, 27%, 39%, and 53% of the demand fullled by solar, wind, and hydro resources under NREL118 S/W, Medium S/W, and High S/W cases, respectively. This outcome suggests that as long as reliability concerns can be mitigated, continuous investment on renewable energy will surely help achieve renewable goals, especially in progressive grids, such as the NREL118. CHAPTER 4. TAMING THE DUCK: CAN STOCHASTIC OPTIMIZATION HELP? 66 0 2000 4000 6000 8000 10000 12000 14000 1/20/24 0:00 1/20/24 1:00 1/20/24 2:00 1/20/24 3:00 1/20/24 4:00 1/20/24 5:00 1/20/24 6:00 1/20/24 7:00 1/20/24 8:00 1/20/24 9:00 1/20/24 10:00 1/20/24 11:00 1/20/24 12:00 1/20/24 13:00 1/20/24 14:00 1/20/24 15:00 1/20/24 16:00 1/20/24 17:00 1/20/24 18:00 1/20/24 19:00 1/20/24 20:00 1/20/24 21:00 1/20/24 22:00 1/20/24 23:00 Megawatts Hydro Natural Gas Solar Wind Other Ren. Curtailment Unmet Demand Over-Generation 0 2000 4000 6000 8000 10000 12000 14000 1/16/24 0:00 1/16/24 1:00 1/16/24 2:00 1/16/24 3:00 1/16/24 4:00 1/16/24 5:00 1/16/24 6:00 1/16/24 7:00 1/16/24 8:00 1/16/24 9:00 1/16/24 10:00 1/16/24 11:00 1/16/24 12:00 1/16/24 13:00 1/16/24 14:00 1/16/24 15:00 1/16/24 16:00 1/16/24 17:00 1/16/24 18:00 1/16/24 19:00 1/16/24 20:00 1/16/24 21:00 1/16/24 22:00 1/16/24 23:00 Megawatts DDD - No reserve Over-Generation (1.8 GW) Renewable Curtailment (6.0 GW) Unmet Demand (3.1 GW) Ren. Curtailment Unmet Demand Over-Generation 0 2000 4000 6000 8000 10000 12000 14000 1/16/24 0:00 1/16/24 1:00 1/16/24 2:00 1/16/24 3:00 1/16/24 4:00 1/16/24 5:00 1/16/24 6:00 1/16/24 7:00 1/16/24 8:00 1/16/24 9:00 1/16/24 10:00 1/16/24 11:00 1/16/24 12:00 1/16/24 13:00 1/16/24 14:00 1/16/24 15:00 1/16/24 16:00 1/16/24 17:00 1/16/24 18:00 1/16/24 19:00 1/16/24 20:00 1/16/24 21:00 1/16/24 22:00 1/16/24 23:00 Megawatts SDS - No Reserve Over-Generation (0.0 GW) Renewable Curtailment (8.6 GW) Unmet Demand (0.6 GW) 0 2000 4000 6000 8000 10000 12000 14000 1/16/24 0:00 1/16/24 1:00 1/16/24 2:00 1/16/24 3:00 1/16/24 4:00 1/16/24 5:00 1/16/24 6:00 1/16/24 7:00 1/16/24 8:00 1/16/24 9:00 1/16/24 10:00 1/16/24 11:00 1/16/24 12:00 1/16/24 13:00 1/16/24 14:00 1/16/24 15:00 1/16/24 16:00 1/16/24 17:00 1/16/24 18:00 1/16/24 19:00 1/16/24 20:00 1/16/24 21:00 1/16/24 22:00 1/16/24 23:00 Megawatts DDD - Medium Reserve Over-Generation (0.0 GW) Renewable Curtailment (7.6 GW) Unmet Demand (1.2 GW) 0 2000 4000 6000 8000 10000 12000 14000 1/16/24 0:00 1/16/24 1:00 1/16/24 2:00 1/16/24 3:00 1/16/24 4:00 1/16/24 5:00 1/16/24 6:00 1/16/24 7:00 1/16/24 8:00 1/16/24 9:00 1/16/24 10:00 1/16/24 11:00 1/16/24 12:00 1/16/24 13:00 1/16/24 14:00 1/16/24 15:00 1/16/24 16:00 1/16/24 17:00 1/16/24 18:00 1/16/24 19:00 1/16/24 20:00 1/16/24 21:00 1/16/24 22:00 1/16/24 23:00 Megawatts SDS - Medium Reserve Over-Generation (0.0 GW) Renewable Curtailment (9.2 GW) Unmet Demand (0.2 GW) Figure 4.8: Over-generation, renewable curtailment, unmet demand, and fullled demand respect to fuel-types (no oil-based generation recorded; based on February 16 data and Medium S/W experiments). 4.4. Notes The quandary of renewable integration, rightfully, worries many power engineers, as solar and wind resources are subject to uncertainty and their outputs are hard-to-predict. This chapter aimed to demonstrate both the positive and negative aspects of the quandary, and evaluated the use of advance decision-making mechanisms to mitigate uncertainty. Our experiments revealed that stochastic modeling and optimization can be a promising choice for achieving reliability in the grids of the future. With further developments, these methodologies could become strong candidates for replacing conventional decision-making tools of the power industry. Chapter 5 Conclusion We described a portfolio of research that addressed dierent challenges in power-system operations planning and stochastic programming. In its essence, our research adopted a multi-faceted approach in strengthening the foun- dations of optimal decision-making under uncertainty, for the prospect of minimum cost and maximum reliability in the power grids of the future. Our UC formulation (Chapter 2) can address problems with seven-times longer lookahead than that of the state-of-the-arts, using much smaller computational eort. This capability can be exploited to handle week-long problems, which would take into account intra-week trends in solar and wind energy, therefore provide decisions with better foresight. The PH-BAB framework (Chapter 3) introduces a new use of the PH algorithm, and leads to a powerful SMIP methodology. This contribution can produce hedged-decisions with reasonable computational eort, and aid the progress towards achieving statistical-optimum in daily operations of the power industry. Finally, Chapter 4 demonstrates the impact of stochastic programming on power-system operations planning, and hopefully, incites enough enthusiasm to encourage similar studies in both the academy and the industry. Before concluding this dissertation, we would like to discuss certain directions of future research that stems from our work. We begin with methodological suggestions pertaining to stochastic programming. (i) For many SCPs, the PH algorithm is able to deliver near-optimal solutions with reasonable eort, but computational performance deteriorates as the algorithm approaches to an optimum. Furthermore, many algorithms take nite amount of time to determine an optimal solution for SLPs (i.e., problems with polyhedral feasible sets and linear objective functions), but the convergence behavior of the PH algorithm for these problems is unclear. Algorithmic modications can be pursued to discover nitely-convergent PH-like algorithms. For instance, it is well-known that Augmented Lagrangian Methods (ALMs) can achieve an optimum for LPs in nite number of iterations. Furthermore, it is easy to see that iterating the rst two steps of the PH algorithm (i.e., optimizing with respect toX and ^ X, respectively) till convergence, leads to an optimal solution for the primal subproblem of ALMs (whereX = ^ X). A hybrid algorithm can bridge the gap between ALMs and the PH algorithm (or, more generally, ADMM), and may lead to optimal solutions in nite time. Moreover, the aforementioned primal iterations can be used for online tuning of the PH algorithm. The behavior of primal and dual iterates can be monitored, and these iterations can be performed to stabilize the sequence of primal iterates. (ii) The generality of the PH-BAB algorithm leads to a wide-variety of applications. First, operations engineers may consider upgrading their mathematical programs to achieve better approximations of systems. For instance, piecewise-linear approximations to generation costs in SUC problems can be replaced with their quadratic coun- terparts. Likewise, convex-approximations to alternating-current power ows can be easily incorporated into SUC problems. Second, statistical or machine learning applications of PH-BAB can be explored. Certain non-convex penalty functions (for inducing sparse decisions) cannot be handled well in many prediction problems. Some of these functions could be modeled with binary variables, which would lead to SMICPs. The PH-BAB algorithm may be a promising contender for addressing such problems. 67 CHAPTER 5. CONCLUSION 68 The nal two suggestions will be related to power-system operations planning. (iii) The value of longer-planning horizons for UC and ED problems should be explored in more detail. Using the experimental framework of Chapter 4, it is now possible to assess the impact of week-long UC problems (with transmission considerations) in much greater depth. Moreover, alternative strategies can be tried to mitigate the computational complexity of the resulting problems. For UC problems, integrality restrictions are only necessary for the day-ahead so that these decisions are implementable. Accordingly, the branch-and-bound process can be stopped, as soon as the commitment schedule for the day-ahead is settled, while the continuous relaxation of the lookahead may still provide important foresight. For both UC and ED problems, varying problem-resolutions can be considered. In particular, the resolutions can be reduced as the planning horizon grows, which would lead to signicant gain in computational performance, and the possibility of considering greater lookahead. More importantly, the solution quality may not be severely aected, as the forecasts for latter periods of a planning horizon are already uncertain, and increased model-complexity for such periods would likely to cause more generalization error. (iv) Finally, the experiments in Chapter 4 should be conducted on dierent test grids, which are more represen- tative of today's power systems. In particular, the considered grid contains no coal or nuclear generators, and the majority of the demand is satised using natural gas resources. This is not very typical. Furthermore, our scenario generation process is independent of the forecasting process, since the latter was simply an input to our studies. This leads low-probability scenarios to be over-represented in stochastic optimization models. Accordingly, the outcomes of our studies might be undermining the value of stochastic optimization in the considered planning prob- lems. Further research will allow us to make more conclusive statements on the impact of stochastic programming on power-system operations. Appendix A Miscellaneous Details on Data Generation for Deterministic UC Problems A.1. Synthetic Instances The instances studied in Ostrowski et al. (2012) consider a time horizon of 24 hours, where the hourly demands are given as percentages of the total generation capacity of the system, and the hourly reserve amounts are set to 3% of the resulting demand. To increase the number of generators, we duplicated all generators tenfold, as in Morales-Espa~ na et al. (2013). To extend the time horizon, we utilized the demand pattern in our realistic instance (see Appendix A.2). For hour h2f1::: 24g and day d 2, we compute d h;d =d h;1 , where d h;d is the corresponding demand in the realistic instance. These ratios are multiplied with the demand percentages given in Ostrowski et al. (2012), and rounded to the nearest integer (see Table A.1). Days Days Hrs 1 2 3 4 5 6 7 Hrs 1 2 3 4 5 6 7 1 71 68 66 64 63 63 63 13 82 87 87 85 83 87 83 2 65 63 60 58 58 57 58 14 80 85 86 84 82 87 82 3 62 60 58 55 55 55 55 15 79 84 86 84 82 87 82 4 60 59 56 54 54 53 52 16 79 84 87 84 82 88 82 5 58 58 55 52 53 52 50 17 83 88 91 88 86 92 85 6 58 60 57 54 56 54 50 18 91 95 97 94 92 97 91 7 60 66 63 59 62 59 52 19 90 93 93 91 89 92 89 8 64 71 68 65 67 64 56 20 88 90 89 88 87 88 86 9 73 80 76 74 75 73 65 21 85 87 85 84 84 83 82 10 80 86 83 81 81 81 75 22 84 85 83 82 82 81 81 11 82 87 86 84 83 84 79 23 79 79 77 76 76 76 77 12 83 88 87 85 84 87 83 24 74 72 71 70 70 70 72 Table A.1: Hourly demands, given as percentages of the total generation capacity of the system (%). A.2. Realistic Instances The realistic instances of §2.3 are based on the winter test problem, made available by FERC 1 and documented in Krall et al. (2012). Here, we only list our assumptions. We set C g to the seasonal capabilities of the generators, 1 http://www.ferc.gov/industries/electric/indus-act/market-planning/rto-commit-test.asp 69 APPENDIX A. DATA GENERATION FOR DETERMINISTIC UC PROBLEMS 70 and let S g = (0:7) R g ; S g = (0:7) R g . C g is either set to the given values, or to minf S g ; S g g, whichever is the minimum. If no positive value is available, we set them to 1. Likewise, missing UT g ; DT g are set to 1. We assumed that the minimum up/downtime requirements are not restrictive at t = 0, and generators incur cold start-up cost after 5 time periods. Generators with no cost entries are neglected (104 cases). Finally, we consider the realized demand data 2 of the PJM region between 01/31/2010 and 01/06/2010, where the former marks the date pertaining to the winter test problem. 2 http://www.pjm.com/markets-and-operations/ops-analysis/historical-load-data.aspx Appendix B Stochastic Server Location Problem SMIP Formulation We provide the formulation of the SSLPs of Ntaimo and Sen (2005), which we use in §3.4. The problem considers a set of locationsJ , zonesZ, and clientsI. The goal is to obtain (i) a cost-ecient allocation of servers into locations inJ , and (ii) a cost-ecient assignment of clients to these servers under random client availability. We list the decision variables of the problem below. z j : 1 if server is located at site j, 0 otherwise y s ij : 1 if client i is served by server at location j under scenario s, 0 otherwise y s j0 : amount of demand over ow due to server capacity limitations. The variable y s j0 must not take large values, therefore is penalized using the penalty function s j . Each problem instance takes the following input parameters: c j : cost of locating a server at location j2J q ij : revenue from client i2I when served by server at location j d ij : demand of client i from server at location j u: server capacity v: upper bound on the total number of located servers w z : minimum number of servers to be located in zone z2Z J z : subset of server locations that belong to zone z i : random availability of client i s i : realization of i under scenario s (1 if client i is present in scenario s, 0 otherwise). The rst- and the second-stage subproblems are respectively given below. min X j2J cjzj +E[h(z; )] h(z; s ) = min X i2I X j2J q s ij y s ij X j2J s j y s j0 s.t. X j2J zjv (B.1a) s.t. X i2I dijy s ij y s j0 uzj; 8i2I; j2J (B.1c) X j2Jz zjwz; 8z2Z (B.1b) X j2J y s ij = s i ; 8i2I (B.1d) zj2B; 8z2Z; y s ij 2B; 8i2I; j2J y s j0 0; 8j2J: Constraints (B.1a) limit the maximum number of server allocations; (B.1b) ensures each zone receives the minimum number of required servers; (B.1c) models the server capacity restrictions; (B.1d) ensures that each client is served by a single server. In the original SSLP instances of Ntaimo and Sen (2005), the penalty function is assumed to have a linear form s j (y s j0 ) =q s j0 y s j0 , where q s j0 is a large number. The quadratic extension of the problem upgrades the penalty 71 APPENDIX B. STOCHASTIC SERVER LOCATION PROBLEM SMIP FORMULATION 72 function as s j (y s j0 ) =q s j0 y s j0 2 . In order to increase the importance of the quadratic penalty function, the server capacity parameter u is adjusted to be smaller in the new instances. Acronyms ADMM alternating direction method of multipliers. ALM Augmented Lagrangian Method. CAISO California ISO. DA-UC day-ahead unit commitment. ED economic dispatch. FERC Federal Energy Regulatory Commission. GHG greenhouse gas. ISO independent system operator. LP linear program, or linear programming. MICP mixed-integer convex program. MILP mixed-integer linear program. MIP mixed-integer program, or mixed-integer pro- gramming. MLR a UC formulation that appeared in Morales- Espa~ na et al. (2013). MWh megawatt hour. NREL National Renewable Energy Lab. OAV a UC formulation that appeared in Ostrowski et al. (2012). PH Progressive Hedging. PH-BAB Progressive-Hedging-based branch-and- bound. PH-BAB-Apx An approximate variant of PH-BAB algorithm. PH-Heuristic The PH algorithm, executed with modied subproblems dened as MICPs (see §3.1). RHA receding (rolling) horizon approach. SCP stochastic convex program. SDM simplicial decomposition method. SED stochastic ED. SLP stochastic linear program. SMICP stochastic mixed-integer convex program. SMILP stochastic mixed-integer linear program. SMIP stochastic mixed-integer program, or stochas- tic mixed-integer programming. SSLP stochastic server location problem. ST-UC short-term unit commitment. STF state-transition formulation. SUC stochastic UC. UC unit commitment. VAR vector autoregression. 73 Bibliography Abrell, J. and Kunz, F. (2015). Integrating Intermittent Renewable Wind Generation - A Stochastic Multi-Market Electricity Model for the European Electricity Market. Networks and Spatial Economics, 15(1):117{147. Ahmed, S. (2013). A scenario decomposition algorithm for 0-1 stochastic programs. Operations Research Letters, 41(6):565{ 569. Ahmed, S., Garcia, R., Kong, N., Ntaimo, L., Parija, G., Qiu, F., and Sen, S. (2015). SIPLIB: A stochastic integer programming test problem library. Angulo, G., Ahmed, S., and Dey, S. S. (2016). Improving the integer L-shaped method. INFORMS Journal on Computing, 28(3):483{490. Arroyo, J. M. and Conejo, A. J. (2000). Optimal response of a thermal unit to an electricity spot market. IEEE Transactions on Power Systems, 15(3):1098{1104. Atakan, S., Lulli, G., and Sen, S. (2018). A state transition MIP formulation for the unit commitment problem. IEEE Transactions on Power Systems, 33(1). Atakan, S. and Sen, S. (2018). A progressive hedging based branch-and-bound algorithm for stochastic mixed-integer programs. Computational Management Science. Preprint. Baldick, R. (1995). The generalized unit commitment problem. IEEE Transactions on Power Systems, 10(1):465{475. Bard, J. F. (1988). Short-term scheduling of thermal-electric generators using Lagrangian relaxation. Operations Research, 36(5):756{766. Barnett, J., Watson, J.-P., and Woodru, D. L. (2017). BBPH: Using progressive hedging within branch and bound to solve multi-stage stochastic mixed integer programs. Operations Research Letters, 45(1):34 { 39. Barth, R., Brand, H., Meibom, P., and Weber, C. (2006). A stochastic unit-commitment model for the evaluation of the impacts of integration of large amounts of intermittent wind power. 2006 9th International Conference on Probabilistic Methods Applied to Power Systems, PMAPS. Benders, J. F. (1962). Partitioning procedures for solving mixed-variables programming problems. Numerische mathematik, 4(1):238{252. Bertsimas, D., Gupta, S., and Lulli, G. (2014). Dynamic resource allocation: A exible and tractable modelling framework. European Journal of Operational Research, 236(1):14{26. Bertsimas, D., Litvinov, E., Sun, X. A., Zhao, J., and Zheng, T. (2013). Adaptive robust optimization for the security constrained unit commitment problem. IEEE Transactions on Power Systems, 28(1):52{63. Birge, J. and Louveaux, F. (1997). Introduction to stochastic programming. Springer, New York. Birge, J. R. (1985). Decomposition and partitioning methods for multistage stochastic linear programs. Operations Research, 33(5):989{1007. Bistline, J. E. and Blanford, G. J. (2016). More than one arrow in the quiver: Why 100% renewables misses the mark. Proceedings of the National Academy of Sciences, 113(28):E3988{E3988. Bixby, R. E. (2002). Solving real-world linear programs: A decade and more of progress. Operations Research, 50(1):3{15. Boland, N., Christiansen, J., Dandurand, B., Eberhard, A., Linderoth, J., Luedtke, J., and Oliveira, F. (2016). Com- bining progressive hedging with a Frank-Wolfe method to compute Lagrangian dual bounds in stochastic mixed-integer programming. Technical report, Optimization Online. Boland, N. L. and Eberhard, A. C. (2014). On the augmented Lagrangian dual for integer programming. Mathematical Programming, 150(2):491{509. Boyd, S. (2010). Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends® in Machine Learning, 3(1):1{122. California ISO (2016). Fast facts: What the duck curve tells us about managing a green grid. Technical report. Retrieved 74 BIBLIOGRAPHY 75 from https://www.caiso.com/Documents/FlexibleResourcesHelpRenewables_FastFacts.pdf. California ISO (2017a). Business practice manual for market operations. Technical report. Retrieved from https://www. caiso.com/rules/Pages/BusinessPracticeManuals/Default.aspx. California ISO (2017b). Fast facts: Impacts of renewable energy on grid operations. Technical report. Retrieved from https://www.caiso.com/Documents/CurtailmentFastFacts.pdf. California ISO (2018a). Flexible resource adequacy criteria and must oer obligations. Retrieved from https://www.caiso. com/informed/Pages/StakeholderProcesses/FlexibleResourceAdequacyCriteria-MustOfferObligations.aspx. California ISO (2018b). Managing oversupply. Retrieved from http://www.caiso.com/informed/Pages/ManagingOversupply. aspx. Care, C. C. and Schultz, R. (1999). Dual decomposition in stochastic integer programming. Operations Research Letters, 24(1-2):37{45. Care, C. C. and Tind, J. (1998). L-shaped decomposition of two-stage stochastic programs with integer recourse. Mathe- matical Programming, 83:451{464. Carpentier, P., Gohen, G., Culioli, J., and Renaud, A. (1996). Stochastic optimization of unit commitment: a new decom- position framework. IEEE Transactions on Power Systems, 11(2):1067{1073. Carri on, M., Arroyo, J. M., and Arroyo, J., M. (2006). A computationally ecient mixed-integer linear formulation for the thermal unit commitment problem. IEEE Transactions on Power Systems, 21(3):1371{1378. Cheung, K., Gade, D., Silva-Monroy, C., Ryan, S. M., Watson, J.-P., Wets, R. J.-B., and Woodru, D. L. (2015). Toward scalable stochastic unit commitment. Energy Systems, pages 1{22. Clack, C. T. M., Qvist, S. A., Apt, J., Bazilian, M., Brandt, A. R., Caldeira, K., Davis, S. J., Diakov, V., Handschy, M. A., Hines, P. D. H., Jaramillo, P., Kammen, D. M., Long, J. C. S., Morgan, M. G., Reed, A., Sivaram, V., Sweeney, J., Tynan, G. R., Victor, D. G., Weyant, J. P., and Whitacre, J. F. (2017). Evaluation of a proposal for reliable low-cost grid power with 100% wind, water, and solar. Proceedings of the National Academy of Sciences, 114(26):6722{6727. Cohen, G. and Zhu, D. L. (1984). Decomposition coordination methods in large scale optimization problems. the nondier- entiable case and the use of augmented Lagrangians. Advances in large scale systems, 1:203{266. Damc-Kurt, P., K u c ukyavuz, S., Rajan, D., and Atamt urk, A. (2015). A polyhedral study of production ramping. Mathe- matical Programming, 158(1):175{205. Denholm, P., Margolis, R., and Milford, J. (2008). Production cost modeling for high levels of photovoltaics penetration. Technical report, National Renewable Energy Laboratory. Retrieved from https://www.nrel.gov/docs/fy08osti/42305. pdf. Denholm, P., O'Connell, M., Brinkman, G., and Jorgenson, J. (2015). Overgeneration from solar energy in California: A eld guide to the duck chart. Technical report, National Renewable Energy Laboratory. Retrieved from http://www.nrel. gov/docs/fy16osti/65023.pdf. Doherty, R. and O'Malley, M. (2005). A new approach to quantify reserve demand in systems with signicant installed wind capacity. IEEE Transactions on Power Systems, 20(2):587{595. Dupa cov a, J., Gr owe-Kuska, N., and R omisch, W. (2003). Scenario reduction in stochastic programming. Mathematical Programming, 95(3):493{511. Escudero, L. F., Gar n, A., Merino, M., and P erez, G. (2009). BFC-MSMIP: An exact branch-and-x coordination approach for solving multistage stochastic mixed 0-1 problems. Top, 17(1):96{122. Eurobarometer (2015). Special eurobarometer 435 \climate change" report. Technical report. Retrieved from https: //ec.europa.eu/clima/sites/clima/files/support/docs/report_2015_en.pdf. European Commission (2009). Renewable energy directive. Retrieved from https://ec.europa.eu/energy/en/topics/ renewable-energy/renewable-energy-directive. Feizollahi, M. J., Ahmed, S., and Sun, A. (2017). Exact augmented Lagrangian duality for mixed integer linear programming. Mathematical Programming, 161:365{387. FERC (2012). Integration of variable energy resources (order no. 764). Federal Energy Regulatory Commission. Retrieved from https://www.ferc.gov/whats-new/comm-meet/2012/062112/E-3.pdf. FERC (2015). Energy primer: A handbook of energy market basics. Technical report, Federal Energy Regulatory Commission. Retrieved from https://www.ferc.gov/market-oversight/guide/energy-primer.pdf. Forsythe, M. (2017). China aims to spend at least $360 billion on renewable energy by 2020. New York Times. Retrieved from https://www.nytimes.com/2017/01/05/world/asia/china-renewable-energy-investment.html?_r=0. Funk, C. and Kennedy, B. (2016). The politics of climate. Technical report, Pew Research Center. Retrieved from http: //assets.pewresearch.org/wp-content/uploads/sites/14/2016/10/10092738/PS_2016.10.04_Politics-of-Climate_ FINAL.pdf. BIBLIOGRAPHY 76 Gade, D., Hackebeil, G., Ryan, S. M., Roger, J.-P. W., and David, J. W. (2016). Obtaining lower bounds from the progressive hedging algorithm for stochastic mixed-integer programs. Mathematical Programming, 157(1):47{67. Gade, D., K u c ukyavuz, S., and Sen, S. (2014). Decomposition algorithms with parametric Gomory cuts for two-stage stochastic integer programs. Mathematical Programming, 144(1-2):39{64. Gangammanavar, H. and Sen, S. (2017). Two-scale stochastic optimization framework for controlling distributed storage devices. IEEE Transactions on Smart Grid. Available as preprint. Gangammanavar, H., Sen, S., and Zavala, V. M. (2015). Stochastic optimization of sub-hourly economic dispatch with wind energy. IEEE Transactions on Power Systems, 31(2):949{959. Garver, L. L. (1962). Power generation scheduling by integer programming - development of theory. Power Apparatus and Systems, Part III, Transactions of the American Institute of Electrical Engineers, 81(3):730{734. Gentile, C., Morales-Espana, G., and Ramos, A. (2016). A tight MIP formulation of the unit commitment problem with start-up and shut-down constraints. EURO Journal on Combinatorial Optimization. Guignard, M. and Kim, S. (1987). Lagrangean Decomposition: A model yielding stronger lagrangean bounds. Mathematical Programming, 39(2):215{228. Guo, G., Hackebeil, G., Ryan, S. M., Watson, J.-P., and Woodru, D. L. (2015). Integration of progressive hedging and dual decomposition in stochastic integer programs. Operations Research Letters, 43(3):311{316. Helmberg, C. (2000). Semidenite programming for combinatorial optimization. Konrad-Zuse-Zentrum f ur Informationstech- nik Berlin. Hering, A. S., Kazor, K., and Kleiber, W. (2015). A markov-switching vector autoregressive stochastic wind generator for multiple spatial and temporal scales. Resources, 4(1):70{92. Higle, J. L. and Sen, S. (1994). Finite master programs in regularized stochastic decomposition. Mathematical Programming, 67(1):143{168. Higle, J. L. and Sen, S. (2006). Multistage stochastic convex programs: Duality and its implications. Annals of Operations Research, 142(1):129{146. Hobbs, B., Rothkopf, M., O'Neill, R., and Chao, H.-P., editors (2001). The Next Generation of Electric Power Unit Commitment Models. International Series in Operations Research & Management Science. Springer-Verlag. Horst, R. and Tuy, H. (2013). Global optimization: Deterministic approaches. Springer Science & Business Media. Jacobson, M. Z., Delucchi, M. A., Cameron, M. A., and Frew, B. A. (2015). Low-cost solution to the grid reliability problem with 100% penetration of intermittent wind, water, and solar for all purposes. Proceedings of the National Academy of Sciences, 112(49):15060{15065. Jiang, R., Guan, Y., and Watson, J. P. (2016). Cutting planes for the multistage stochastic unit commitment problem. Mathematical Programming, 157(1):121{151. J ornsten, K., N asberg, M., and Smeds, P. (1985). Variable splitting: A new Lagrangean relaxation approach to some mathematical programming models. Technical report, University of Link oping, Department of Mathematics. Kall, P. and Wallace, S. W. (1994). Stochastic Programming. John Wiley & Sons. Krall, E., Higgins, M., and O'Neill, R. P. (2012). RTO unit commitment test system. Technical report, Federal Energy Regulatory Commission. Retrieved from http://www.ferc.gov/legal/staff-reports/rto-COMMITMENT-TEST.pdf. Laporte, G. and Louveaux, F. V. (1993). The integer L-shaped method for stochastic integer programs with complete recourse. Operations Research Letters, 13(3):133{142. Lee, J., Leung, J., and Margot, F. (2004). Min-up / min-down polytopes. Discrete Optimization, 1:77{85. Lee, Y. Y. and Baldick, R. (2013). A frequency-constrained stochastic economic dispatch model. IEEE Transactions on Power Systems, 28(3):2301{2312. Li, T. and Shahidehpour, M. (2005). Price-based unit commitment: a case of Lagrangian relaxation versus mixed integer programming. IEEE Transactions on Power Systems, 20(4):2015{2025. Lkketangen, A. and Woodru, D. L. (1996). Progressive hedging and tabu search applied to mixed integer (0,1) multistage stochastic programming. Journal of Heuristics, 2(2):111{128. Louveaux, F. V. (1980). A solution method for multistage stochastic programs with recourse with application to an energy investment problem. Operations Research, 28(4):889{902. Lubin, M., Martin, K., Petra, C. G., and Sandk c, B. (2013). On parallelizing dual decomposition in stochastic integer programming. Operations Research Letters, 41(3):252 { 258. Lulli, G. and Sen, S. (2004). A branch-and-price algorithm for multistage stochastic integer programming with application to stochastic batch-sizing problems. Management Science, 50(6):786{796. Lutkepohl, H. (2007). New Introduction to Multiple Time Series Analysis. Springer Publishing Company Inc. Makarov, Y. V., Huang, Z., Etingov, P. V., Ma, J., Guttromson, R. T., Subbarao, K., and Chakrabarti, B. B. (2010). BIBLIOGRAPHY 77 Incorporating wind generation and load forecast uncertainties into power grid operations. Technical Report PNNL-19189, Pacic Northwest National Laboratory. Monteiro, C., Bessa, R., Botterud, A., Wang, J., and Conzelmann, G. (2009). Wind power forecasting: State-of-the-art 2009. Technical report, Argonne National Laboratory. Retrieved from https://www.anl.gov/energy-systems/publication/ wind-power-forecasting-state-art-2009. Morales, J., M nguez, R., and Conejo, A. J. (2010). A methodology to generate statistically dependent wind speed scenarios. Applied Energy, 87(3):843 { 855. Morales, J. M., Conejo, A. J., and Perez-Ruiz, J. (2009). Economic valuation of reserves in power systems with high penetration of wind power. IEEE Transactions on Power Systems, 24(2):900{910. Morales-Espa~ na, G., Latorre, J. M., and Ramos, A. (2013). Tight and compact MILP formulation for the thermal unit commitment problem. IEEE Transactions on Power Systems, 28(4):4897{4908. Muckstadt, J. and Koenig, S. A. (1977). An application of Lagrangian relaxation to scheduling in power-generation systems. Operations Research, 25(3):387{403. Mulvey, J. M. and Vladimirou, H. (1991). Applying the progressive hedging algorithm to stochastic generalized networks. Annals of Operations Research, 31(1):399{424. Mulvey, J. R. and Ruszczy nski, A. (1995). A new scenario decomposition method for large-scale stochastic optimization. Operations Research, 43(3):477{490. Nemhauser, G. L. and Wolsey, L. A. (1999). Integer programming and combinatorial optimization. John Wiley and Sons, New York. New York ISO (2016). Power trends 2016: The changing energy landscape. Technical report. Retrieved from http://www.nyiso.com/public/webdocs/media_room/publications_presentations/Power_Trends/Power_Trends/ 2016-power-trends-FINAL-070516.pdf. New York ISO (2017a). Day-ahead scheduling manual. Technical report. Retrieved from http://www.nyiso.com/public/ webdocs/markets_operations/documents/Manuals_and_Guides/Manuals/Operations/dayahd_schd_mnl.pdf. New York ISO (2017b). Power trends 2017: New York's evolving electric grid. Technical report. Retrieved from http: //www.nyiso.com/public/webdocs/media_room/publications_presentations/Power_Trends/Power_Trends/2017_Power_ Trends.pdf. Nikolewski, R. (2017). Can California really hit a 100% renewable energy target? San Diego Tribune. Retrieved from http: //www.sandiegouniontribune.com/business/energy-green/sd-fi-california-100percent-20170601-story.html. Nowak, M. P. and R omisch, W. (2000). Stochastic Lagrangian relaxation applied to power scheduling in a hydro-thermal system under uncertainty. Annals of Operations Research, 100(1-4):251{272. Ntaimo, L. and Sen, S. (2005). The million-variable \march" for stochastic combinatorial optimization. Journal of Global Optimization, 32(3):385{400. Obama, B. (2017). The irreversible momentum of clean energy. Science. Obama, B. H. (2013). Federal Leadership on Energy Management [Presidential Memorandum]. Washington, DC: The White House. Retrieved from https://obamawhitehouse.archives.gov/the-press-office/2013/12/05/ presidential-memorandum-federal-leadership-energy-management. Olson, A., Mahone, A., Hart, E., Hargreaves, J., Jones, R., Schlag, N., Kwok, G., Ryan, N., Orans, R., and Frowd, R. (2015). Halfway There: Can California Achieve a 50% Renewable Grid? IEEE Power and Energy Magazine, 13(4):41{52. O'Neill, R. P., Dautell, T., and Krall, E. (2011). Recent ISO software enhancements and future software and modeling plans. Technical report, Federal Energy Regulatory Commission. Retrieved from http://www.ferc.gov/industries/electric/ indus-act/rto/rto-iso-soft-2011.pdf. Ostrowski, J., Anjos, M. F., and Vannelli, A. (2012). Tight mixed integer linear programming formulations for the unit commitment problem. IEEE Transactions on Power Systems, 27(1):39{46. Ozturk, U., Mazumdar, M., and Norman, B. (2004). A solution to the stochastic unit commitment problem using chance constrained programming. IEEE Transactions on Power Systems, 19(3):1589{1598. Padhy, N. (2004). Unit commitment - A bibliographical survey. IEEE Transactions on Power Systems, 19(2):1196{1205. Pages-Bernausa, A., P erez-Vald es, G., and Tomasgard, A. (2015). A parallelised distributed implementation of a branch and x coordination algorithm. European journal of operational research, 244:77{85. Pan, K. and Guan, Y. (2016). A polyhedral study of the integrated minimum-up/-down time and ramping polytope. Optimization Online. Pang, C. K., Sheble, G. B., and Albuyeh, F. (1981). Evaluation of dynamic programming based methods and multiple area representation for thermal unit commitments. IEEE Transactions on Power Apparatus and Systems, 100(3):1212{1218. Papavasiliou, A. and Oren, S. S. (2013). Multiarea stochastic unit commitment for high wind penetration in a transmission BIBLIOGRAPHY 78 constrained network. Operations Research, 61(3):578{592. Papavasiliou, A., Oren, S. S., and Neill, R. P. O. (2011). Reserve requirements for wind power integration: A stochastic programming framework. IEEE Transactions on Power Systems, 26(4):2197{2206. Pena, I., Brancucci, C., and Hodge, B. M. (2017). An Extended IEEE 118-bus Test System with High Renewable Penetration. IEEE Transactions on Power Systems, 33(1):281{289. Penn, I. (2017). California invested heavily in solar power. now there's so much that other states are sometimes paid to take it. Los Angeles Times. Retrieved from http://www.latimes.com/projects/la-fi-electricity-solar/. Pritsker, A. A. B., Waiters, L. J., and Wolfe, P. M. (1969). Multiproject scheduling with limited resources: A zero-one programming approach. Management Science, 16(1):93{108. Qi, Y. and Sen, S. (2016). The Ancestral Benders' cutting plane algorithm with multi-term disjunctions for mixed-integer recourse decisions in stochastic programming. Mathematical Programming. Rajan, D. and Takriti, S. (2005). Minimum up / down polytopes of the unit commitment problem with start-up costs. Technical report, IBM. Ralphs, T. and Hassanzadeh, A. (2014). A generalization of Benders' algorithm for two-stage stochastic optimization problems with mixed integer recourse. Technical report, Laboratory for Computational Optimization at Lehigh (COR@L), Department of Industrial and Systems Engineering, Lehigh University. Rios, I., Wets, R. J.-B., and Woodru, D. L. (2015). Multi-period forecasting and scenario generation with limited data. Computational Management Science, 12(2):267{295. Rockafellar, R. T. (1976). Monotone operators and the proximal point algorithm. SIAM Journal on Control and Optimization, 14(5):877{898. Rockafellar, R. T. and Wets, R. J.-B. (1991). Scenarios and policy aggregation in optimization under uncertainty. Mathematics of Operations Research, 16(1):119{147. Ruiz, P. A., Philbrick, C. R., and Sauer, P. W. (2009). Wind power day-ahead uncertainty management through stochastic unit commitment policies. 2009 IEEE/PES Power Systems Conference and Exposition, PSCE 2009, pages 1{9. Ruszczy nski (2002). Probabilistic programming with discrete distributions and precedence constrained knapsack polyhedra. Mathematical Programming, 93:195{215. Ryan, K., Ahmed, S., Dey, S. S., and Rajan, D. (2016). Optimization driven scenario grouping. Technical report, School of Industrial & Systems Engineering, Georgia Institute of Technology. Ryan, S. M., Wets, R. J. B., Woodru, D. L., Silva-Monroy, C., and Watson, J. P. (2013). Toward scalable, parallel progressive hedging for stochastic unit commitment. In IEEE Power Energy Society (PES) General Meeting, pages 1{5. Saravanan, B., Das, S., Sikri, S., and Kothari, D. P. (2013). A solution to the unit commitment problem - a review. Frontiers in Energy, 7:223{236. SB-100 (2018). California Renewables Portfolio Standard Program: emissions of greenhouse gases. California State Bill, (in assembly). SB-350 (2015). Clean Energy and Pollution Reduction Act of 2015. California State Bill, Chapter 547. Sen, S. (2005). Algorithms for Stochastic Mixed-Integer Programming Models. In Aardal, K., Nemhauser, G. L., and Weismantel, R., editors, Handbooks in Operations Research and Management Science, chapter 9, pages 515{558. Elsevier. Sen, S. and Deng, Y. (2018). Learning enabled optimization: Towards a fusion of statistical learning and stochastic pro- gramming. Technical report, University of Southern California. Available in http://www.optimization-online.org/DB_ FILE/2017/03/5904.pdf. Sen, S. and Higle, J. L. (2005). The C 3 theorem and a D 2 algorithm for large scale stochastic mixed-integer programming: Set convexication. Mathematical Programming, 104(1):1{20. Sen, S. and Liu, Y. (2016). Mitigating uncertainty via compromise decisions in two-stage stochastic linear programming: Variance reduction. Operations Research, 64(6):1422{1437. Sen, S. and Sherali, H. D. (2006). Decomposition with branch-and-cut approaches for two-stage stochastic mixed-integer programming. Mathematical Programming, Series A, 106(2):203{223. Sen, S., Yu, L., and Genc, T. (2006). A stochastic programming approach to power portfolio optimization. Operations Research, 54(1):55{72. Shahidehpour, M., Yamin, H., and Li, Z. (2002). Market Operations in Electric Power Systems: Forecasting, Scheduling, and Risk Management. IEEE, Wiley-Interscience. Shiina, T. and Birge, J. R. (2004). Stochastic unit commitment problem. International Transactions in Operational Research, 11(1):19{32. Takriti, S. and Birge, J. (2000). Using integer programming to rene Lagrangian-based unit commitment solutions. IEEE Transactions on Power Systems, 15(1):151{156. BIBLIOGRAPHY 79 Takriti, S., Birge, J. R., and Long, E. (1996). A stochastic model for the unit commitment problem. IEEE Transactions on Power Systems, 11(3):1497{1508. The Economist (2017). Renewable energy: a world upside down. Retrieved from http://www.economist.com/news/briefing/ 21717365-wind-and-solar-energy-are-disrupting-century-old-model-providing-electricity-what-will. The Guardian (2017). 'spectacular' drop in renewable energy costs leads to record global boost. Retrieved from https : / / www . theguardian . com / environment / 2017 / jun / 06 / spectacular-drop-in-renewable-energy-costs-leads-to-record-global-boost. Trapp, A. C., Prokopyev, O. A., and Schaefer, A. J. (2013). On a level-set characterization of the value function of an integer program and its application to stochastic programming. Operations Research, 61(2):498{511. Tseng, C., Li, C., and Oren, S. (2000). Solving the unit commitment problem by a unit decommitment method. Journal of Optimization Theory, 105:707{730. Tuohy, A., Zack, J., Haupt, S. E., Sharp, J., Ahlstrom, M., Dise, S., Grimit, E., Mohrlen, C., Lange, M., Casado, M. G., Black, J., Marquis, M., and Collier, C. (2015). Solar forecasting: Methods, challenges, and performance. IEEE Power and Energy Magazine, 13(6):50{59. Turgeon, A. (1977). Optimal unit commitment. IEEE Transactions on Automatic Control, 23(2):223{227. Van Slyke, R. M. and Wets, R. (1969). L-shaped linear programs with applications to optimal control and stochastic programming. SIAM Journal on Applied Mathematics, 17(4):638{663. Wang, J., Shahidehpour, M., and Li, Z. (2008). Security-constrained unit commitment with volatile wind power generation. IEEE Transactions on Power Systems, 23(3):1319{1327. Watson, J.-P. and Woodru, D. L. (2010). Progressive hedging innovations for a class of stochastic mixed-integer resource allocation problems. Computational Management Science, 8(4):355{370. Wollmer, R. D. (1980). Two stage linear programming under uncertainty with 0-1 integer rst stage variables. Mathematical Programming, 19(1):279{288. Woodru, D. L., Deride, J., Staid, A., Watson, J.-P., Slevogt, G., and Silva-Monroy, C. (2018). Constructing probabilistic scenarios for wide-area solar power generation. Solar Energy, 160:153 { 167. Zavala, V. M., Kim, K., Anitescu, M., and Birge, J. (2017). A stochastic electricity market clearing formulation with consistent pricing properties. Operations Research, 65(3):557{576. Zhang, M. and K u c ukyavuz, S. (2014). Finitely convergent decomposition algorithms for two-stage stochastic pure integer programs. SIAM Journal on Optimization, 24(4):1933{1951. Zheng, Q. P., Wang, J., Pardalos, P. M., and Guan, Y. (2013). A decomposition approach to the two-stage stochastic unit commitment problem. Annals of Operations Research, 210(1):387{410.

Abstract (if available)

Abstract The momentum of renewable energy is transforming the power grids across the world. A markedly increasing share of electricity-demand is sourced by solar and wind resources, which are volatile and difficult-to-store, therefore causing concern regarding the reliability of power-system operations. This thesis adopts a multi-faceted approach on expanding our knowledge of decision-making tools which could eventually be beneficial to mitigate this uncertainty, through improved optimization of day-to-day operations. ❧ We begin with the so-called unit-commitment problem and propose a novel formulation that is tailored for problems with longer planning horizons. With its strong polyhedral properties, our formulation can deliver optimal solutions to week-long and industry-strength problems in nearly 10 minutes, using an ordinary desktop computer. Furthermore, our research shows that considering week-long planning horizons, in lieu of the industry-standard day-long ones, captures demand trends better, and promises improved operating-costs. Likewise, this formulation can accommodate the intra-week trends in solar and wind resources, and may produce decisions with better foresight. ❧ Our second contribution is a methodological one that aims to address uncertainty in general mixed-integer formulations. In particular, we consider the well-known Progressive Hedging (PH) algorithm from the stochastic convex optimization literature, and extend it for solving multi-stage stochastic mixed-integer convex programs (SMICPs). Our framework for using PH is guaranteed to compute optimal solutions for SMICPs, which contrasts with most previous extensions of PH, as they have been implemented without convergence guarantees. More importantly, we demonstrate the effectiveness of our framework through computational experiments, and observe significant improvements over the performances of benchmark algorithms. In particular, our algorithms can deliver optimal solutions to certain instances in nearly 30 seconds, whereas a decade ago, multiple hours of computational time had not been enough to identify an optimal decision. ❧ Finally, we turn our attention to renewable energy, and assess the risks of increasing levels of renewable-integration on power grids. To this end, we develop an experimental framework that resembles the daily operations of power-management authorities. We also use this framework to assess the use of deterministic and stochastic optimization methodologies for different components of power-system planning processes. Our experiments provide quantitative support to the notion that introducing significant renewable resources into power systems requires a systems-view of supply and demand uncertainty, the transmission capacity of the grid, and just as importantly, the operations planning strategy used. With regards to the latter point, we demonstrate that stochastic programming can be a valuable decision-making tool for maintaining the reliability and cost-effectiveness of the grid.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Computational stochastic programming with stochastic decomposition

PDF

Learning and control in decentralized stochastic systems

PDF

Computational validation of stochastic programming models and applications

PDF

I. Asynchronous optimization over weakly coupled renewal systems

PDF

On the interplay between stochastic programming, non-parametric statistics, and nonconvex optimization

PDF

The fusion of predictive and prescriptive analytics via stochastic programming

PDF

Learning enabled optimization for data and decision sciences

PDF

Smarter markets for a smarter grid: pricing randomness, flexibility and risk

PDF

Models and algorithms for the freight movement problem in drayage operations

PDF

The smart grid network: pricing, markets and incentives

PDF

Variation-aware circuit and chip level power optimization in digital VLSI systems

PDF

Empirical methods in control and optimization

PDF

A stochastic conjugate subgradient framework for large-scale stochastic optimization problems

PDF

Mixed-integer nonlinear programming with binary variables

PDF

Landscape analysis and algorithms for large scale non-convex optimization

PDF

Online learning algorithms for network optimization with unknown variables

PDF

Topics in algorithms for new classes of non-cooperative games

PDF

Stochastic games with expected-value constraints

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Iterative path integral stochastic optimal control: theory and applications to motor control

Asset Metadata

Creator Atakan, Semih (author)

Core Title The next generation of power-system operations: modeling and optimization innovations to mitigate renewable uncertainty

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Industrial and Systems Engineering

Publication Date 08/09/2020

Defense Date 07/08/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag economic dispatch,mixed-integer programming,OAI-PMH Harvest,optimization,progressive hedging,renewable energy,stochastic mixed-integer programming,unit commitment

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Sen, Suvrajeet (committee chair), Higle, Julia Lynne (committee member), Jain, Rahul (committee member)

Creator Email atakan@usc.edu,semihatakan@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-67685

Unique identifier UC11671751

Identifier etd-AtakanSemi-6722.pdf (filename),usctheses-c89-67685 (legacy record id)

Legacy Identifier etd-AtakanSemi-6722.pdf

Dmrecord 67685

Document Type Dissertation

Format application/pdf (imt)

Rights Atakan, Semih

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA