Applied Science, Faculty of
Civil Engineering, Department of
DSpace
UBCV
Mitchell, Kyle Terrence
Master of Applied Science - MASc
University of British Columbia
Stochastic Dual Dynamic Programming is a heuristic stochastic optimization algorithm used to analyze multistage decision making problems with stochastic variables through Benders decomposition and solvable linear approximations. Benders cuts are utilized to develop a piecewise linear approximation of a 'benefit-to-go' or future value function at each stage. This thesis assesses BC Hydro's Columbia and Peace River hydroelectric network representation in a multistage stochastic optimization problem with a finite planning horizon. The model attempts to maximize the 'water value' across the MCA, GMS, and ARD facilities over a duration of five years. Existing model parameter-space is heavily restricted due to the use of static market prices and load demands. Changes have been proposed to expand on the current parameter-space pertaining to inflow dependent market electricity prices through the implementation of seasonally forecasted price and load blocks. Preliminary results indicate the existing model implementation does not provide sufficient parameter-space coverage with respect to market prices, resulting in possible over-estimations of water value as the model approaches convergence.
PARAMETER-SPACE DEPENDENCE IN A STOCHASTIC DUAL DYNAMICPROGRAMMINGMODEL OF A MULTIRESERVOIR HYDROPOWER SYSTEMbyKyle Terrence MitchellA THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FORTHE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Civil Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)January 2023© Kyle Terrence Mitchell, 2023iiThe following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the thesis entitled:Parameter-Space Dependence in a Stochastic Dual Dynamic Programming Model of a Multireser-voir Hydropower Systemsubmitted by Kyle Mitchell in partial fulfilment of the requirements forthe degree of Master of Applied Sciencein Civil EngineeringExamining Committee:Ziad Shawwash, Assistant Professor, Civil Engineering, UBCSupervisorGregory Lawrence, Professor, Civil Engineering, UBCSupervisory Committee MemberiiiAbstractStochastic Dual Dynamic Programming is a heuristic stochastic optimization algorithm used toanalyze multistage decision making problems with stochastic variables through Benders decom-position and solvable linear approximations. Benders cuts are utilized to develop a piecewiselinear approximation of a ’benefit-to-go’ or future value function at each stage. This thesis as-sesses BCHydro’s Columbia and Peace River hydroelectric network representation in a multistagestochastic optimization problem with a finite planning horizon. The model attempts to maximizethe ’water value’ across the MCA, GMS, and ARD facilities over a duration of five years. Exist-ing model parameter-space is heavily restricted due to the use of static market prices and loaddemands. Changes have been proposed to expand on the current parameter-space pertaining toinflow dependent market electricity prices through the implementation of seasonally forecastedprice and load blocks. Preliminary results indicate the existing model implementation does notprovide sufficient parameter-space coverage with respect to market prices, resulting in possibleover-estimations of water value as the model approaches convergence.ivLay SummaryStochasticDualDynamicProgramming (SDDP) is apowerfulmathematical tool that helpsdecision-makers solve complex problems that involve multiple stages and uncertainties about the future.It is a type of dynamic programming algorithm, which is a method for breaking down a complexproblem into smaller, more manageable sub-problems. The SDDP algorithm uses both determin-istic (non-random) and stochastic (random) information to make decisions, which allows it to takeinto account the uncertainty of future events. This is particularly useful for problems with a largenumber of decision stages, such as those found in fields like finance, energy, and transportation.SDDP can be used to find the best decision at each stage of the problem, taking into account theuncertainty of future events, and ultimately providing an optimal solution. This thesis appliesSDDP to a real-world hydroelectric system consisting of several interconnected dams.vPrefaceThis thesis is heavily based on an internally developed application of Stochastic Dual DynamicProgramming by BC Hydro. None of the text in this thesis is taken directly from previouslypublished or collaborative articles.The original (existing) model referenced frequently in Chapters 3, 4, and 5 was constructed overseveral years by Z. Shawwash in collaboration with other doctorate and post-doctorate studentsfrom BC Hydro. Additional collaborators have also continuously contributed to this model devel-opment.Any code developed within this thesis is my own work, however, it builds off of concepts intro-duced by Z. Shawwash, specifically with regard to price and load block shaping. Thus, my workconducted in Chapters 3 and 4 should be considered an extension of that originally developed byZ. Shawwash and others involved from BC Hydro.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Importance of Water Resource Use . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Challenges Addressed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Background and Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Theory Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Essential Understandings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.1 Feasibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.2 Boundedness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.3 Optimality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8vii2.2.4 Duality Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3 Decomposable Structures and Benders Decomposition . . . . . . . . . . . . . . . . . 122.4 Optimization using Benders Decomposition . . . . . . . . . . . . . . . . . . . . . . . . 152.4.1 Dual Dynamic Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4.2 Stochastic Dual Dynamic Programming . . . . . . . . . . . . . . . . . . . . . . 202.4.3 Multi-stage Stochastic Problems and Nested Benders Decomposition . . . . . 222.4.4 Convergence Proof and Criterion of SDDP Models . . . . . . . . . . . . . . . . 262.5 Applications from Literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.5.1 Optimal Scheduling of a Hydrothermal Generating System by Pereira & Pinto 272.5.2 Optimal Operation in Deregulated Electricity Markets by Iliadis et al. . . . . . 303 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.1 Model of the BC Hydropower System . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.1.1 Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34Inflow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36Flood Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Treaty Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Non-Treaty Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40Water Mass Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.2 Model Implementation in Julia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2.1 SDDP Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2.2 Model Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.2.3 Price and Load Blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47SDDP Price Blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50SDDP Load Blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.2.4 Price and Load Block Implementation in the Julia Language . . . . . . . . . . 52Existing SDDP Julia Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52Computational Tractability and the Master Sheet . . . . . . . . . . . . . . . . . 523.3 Existing Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.3.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54viii3.3.2 Storage Accounts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.3.3 Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.3.4 Implications of 27 vs 49 Simulation Years . . . . . . . . . . . . . . . . . . . . . 653.3.5 Computational Effort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.1 Implemented Changes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.1.1 Challenges Encountered . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.2 Objective Function Trends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.2.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.3 Storage Accounts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.4 Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.5 Validation of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.5.1 Storage Accounts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.5.2 Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.6 Computational Effort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905 Conclusions and Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.1 Discussion of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.2 Recommendations for Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94ixList of Tables3.1 1972 GOM Load Data (MW) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.2 1972 GOM Price Data (US$/MWh) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.3 1972 Price Blocks (US$/MWh) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.4 1972 Load Blocks (MWh) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5 ε Variations for 49 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.6 ε by Month for 49 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.7 GMS and MCAMarginal Value Trends . . . . . . . . . . . . . . . . . . . . . . . . . . . 64xList of Figures1.1 Top: Projected Share of New Vehicle Sales. Bottom: Projected Increase in Electricityuse in the Transportation Sector (C.E.R., 2018) . . . . . . . . . . . . . . . . . . . . . . . 31.2 Non-convex Function with Benders Cut (Fullner & Rebennack, 2021) . . . . . . . . . 42.1 A Feasible Region Enclosed by Three Constraints (Murphy, 2013) . . . . . . . . . . . 82.2 FromLeft to Right: a Bounded Finite Feasible Space, a Bounded Semi-Infinite Space,an Unbounded Semi-Infinite Space (Murphy, 2013) . . . . . . . . . . . . . . . . . . . 92.3 Optimality of a Linear System in a Finite Bounded Space (Murphy, 2013) . . . . . . . 92.4 Algorithm for Generating Possible πi Values for π∗ . . . . . . . . . . . . . . . . . . . 202.5 Three-stage Stochastic Scenario Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.6 A Typical Stochastic Decision Tree (Murphy, 2013) . . . . . . . . . . . . . . . . . . . . 232.7 Nested Two-stage Problems (Kazempour, 2022) . . . . . . . . . . . . . . . . . . . . . . 232.8 A Three-stage Stochastic Problem (Murphy, 2013) . . . . . . . . . . . . . . . . . . . . 242.9 Expected Operations Cost (Pereira & Pinto, 1991) . . . . . . . . . . . . . . . . . . . . . 302.10 Optimal Point of FBF and IBF Sum (Iliadis et al., 2008) . . . . . . . . . . . . . . . . . . 312.11 Optimal (12 Month) Hydropower Generation Results (Iliadis et al., 2008) . . . . . . . 323.1 Network of Primary BC Hydro Hydroelectric Facilities (Guan et al., 2018) . . . . . . 343.2 Relationships Between Implemented Models (Guan et al., 2018) . . . . . . . . . . . . 353.3 Use of Statistical Methods in Inflow RV Development (Guan et al., 2018) . . . . . . . 373.4 Modeled Reservoirs and Storage Accounts (Guan et al., 2018) . . . . . . . . . . . . . . 393.5 Benders Cuts and water value Approximation . . . . . . . . . . . . . . . . . . . . . . 443.6 Price and Load Blocks Flowchart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.7 Current Price and Load Block Implementation with Model Flow . . . . . . . . . . . . 53xi3.8 Final Simulation Objective Value Trends . . . . . . . . . . . . . . . . . . . . . . . . . . 553.9 ARD Storage Trends (49 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.10 GMS Storage Trends (49 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.11 MCA Storage Trends (49 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.12 ARD Treaty Storage Trends (49 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . 603.13 MCA Treaty Storage Trends (49 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . 603.14 1972 Price Block . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.15 1972 Load Block . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.16 ARD Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.17 GMS Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.18 MCAMarginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.19 GMS Storage Trends (27 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.20 ARD Storage Trends (27 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.21 MCA Storage Trends (27 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.22 ARD Treaty Storage Trends (27 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . . 673.23 MCA Treaty Storage Trends (27 Sim Years) . . . . . . . . . . . . . . . . . . . . . . . . 673.24 Computation Time for Iteration Requirements . . . . . . . . . . . . . . . . . . . . . . 683.25 Computational Time Projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.1 Model After Implemented Changes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.2 New Price Parameter-Space vs. 1972 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.3 New Load Parameter-Space vs. 1972 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.4 Negative Load in LLH Blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.5 New Model Final Simulation Objective Value Trends . . . . . . . . . . . . . . . . . . . 754.6 New Model GMS Storage Trends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.7 New Model ARD Storage Trends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.8 New Model MCA Storage Trends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.9 New Model ARD Treaty Storage Trends . . . . . . . . . . . . . . . . . . . . . . . . . . 784.10 New Model MCA Treaty Storage Trends . . . . . . . . . . . . . . . . . . . . . . . . . . 794.11 New Model GMS Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . 80xii4.12 New Model ARD Marginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . 804.13 New Model MCAMarginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . . . 814.14 New Model TY ARDMarginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . 814.15 New Model TY MCAMarginal Water Value . . . . . . . . . . . . . . . . . . . . . . . . 824.16 SLP vs New SDDP GMS Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.17 SLP vs New SDDP ARD Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.18 SLP vs New SDDP MCA Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.19 SLP vs New SDDP TY ARD Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.20 SLP vs New SDDP TY MCA Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.21 SLP vs New SDDP GMS Marginal Value . . . . . . . . . . . . . . . . . . . . . . . . . . 874.22 SLP vs New SDDP ARDMarginal Value . . . . . . . . . . . . . . . . . . . . . . . . . . 884.23 SLP vs New SDDP MCAMarginal Value . . . . . . . . . . . . . . . . . . . . . . . . . 884.24 SLP vs New SDDP TY ARDMarginal Value . . . . . . . . . . . . . . . . . . . . . . . . 894.25 SLP vs New SDDP TY MCAMarginal Value . . . . . . . . . . . . . . . . . . . . . . . 894.26 New and Old Computation Time by Iteration . . . . . . . . . . . . . . . . . . . . . . . 91xiiiList of AbbreviationsARD Hugh Keenleyside DamBCH British Columbia HydroCRT Columbia River TreatyFBF Future Benefit FunctionGMS GordonMerritt Shrum DamGSO Generation System OperationsIBF Immediate Benefit FunctionLCA Libby Coordination AgreementLP Linear ProgrammingMCA MiCA DamMCS Monte Carlo SamplingPCA Principal Component AnalysisPCN Peace CanyoN DamREV REVelstoke DamRV Random Variable(s)SDDP Stochastic Dual Dynamic ProgrammingSLP Successive Linear ProgrammingSOA Supplementary Operating AgreementSRD Storage Reservation DiagramsTVF Terminal Value FunctionWVM Water ValueModelxivAcknowledgementsThis thesis would not have been possible without the incredible knowledge, support, and in-sight provided by Dr. Shawwash and the corresponding team at BC Hydro lead by Tim Blair.Thank you to the entire team for always providing a helping hand when I was stuck on a conceptor idea.The funding and academic resources provided by UBC were essential in facilitating my continuedenrollment during the Covid-19 pandemic. The quick transition from in-person to remote learningensured students such asmyself were able to graduate on timewithout any significant interruptionor reduction in course quality. So, to the staff and professors who worked tirelessly to make thishappen, thank you.As well, to Dr. Yonemitsu and Dr. Isaacson, I owe a great deal of gratitude for providing excellentTeaching Assistant and learning opportunities throughout my time at UBC. While not directly re-lated to the content of this thesis, the experience and passion I’ve subsequently gained for helpingother students succeed has been invaluable.To my immediate family and partner I owe a great deal. Their patience, love, and support haskept me determined even on the worst of days. Thank you to Miki ミキ and Tae たえ for sharingtheir strength through the past two years. As well, this would not have been possible withoutthe unyielding support of my mother Ana and her partner Simon. Their shared love of sciencesparked my own and undoubtedly contributed to where I am today.1Chapter 1Introduction and MotivationSystem optimization, whether through algorithms or heuristics, is an incredibly broad field ofstudy with new techniques and strategies continuously under study and development around theworld. It is easy to consider these techniques a consequence of the power of computing; however,optimization has been a driving force behind scientific discovery since the conception of geometryby Euclid of Alexandria in 300BC (Euclid, 1996). Whether optimization is used to calculateminimum distance or to determine a shortest path with the modern A* algorithm ((Brewka, 1996),the sum total of these techniques has helped shape virtually every aspect of modern science. Whileone of the earliest models of dynamic programming was developed by Johannes Kepler in 1615to assist in finding himself a new wife (Ferguson, 1989), the application of modern day DynamicProgramming discussed within this thesis pertains to a more resourceful problem, specificallythe behaviours of a Stochastic Dual Dynamic Programming (SDDP) model of a multireservoirhydropower system. This particular SDDP model attempts to maximize the net benefit acrossseveral of BC Hydro’s hydroelectric facilities by generating a water value function dependent onseasonal forecasts and storage policies. The primary focus of this thesis is to test the dependenceof this model on two of its key parameters, namely energy prices and load demands.1.1 Importance of Water Resource UseOver the past several decades, water consumption has more than doubled and imposed increasingglobal water stress levels (Wada et al., 2011). A summary of 10 academic studies from 1999 to2011 by Wada et al., 2011 showed that, when assessing global water stress, anywhere between 8%2and 46% of the global population could be under high water stress at any given time. High waterstress is defined by aWater Stress Index (WSI) of greater than 0.4 (ratio of freshwater availability todemand). From 1960-2005, global levels of high water stress surged from 9% to 35% (Wada et al.,2011).Water is critical to every day life, yet is vastly taken for granted by Western populations. Fromirrigation for crops to generating electricity, a large proportion of water consumption occurs outof view and can be easy to overlook. Additional global trends towards renewable energy shedslight on conscious water consumption as a source of consistent, ’storable’ green energy. RenewableEnergy Sources (RES) within highly developed countries such as Canada have significant potentialto contribute to their economic, social, and environmental sustainability. Canada is one of theleading producers of hydroelectric energy, producing 26% of the national demand in 2019 (Firsovaet al., 2019). Current global leaders in hydroelectric dependence are Norway (64%), Brazil (35%),and Sweden (30%). An analysis of global energy trends by Firsova et al., 2019 concluded thatRES are the strategic future of humanity due to the eventual exhaustion of natural resources,including the long term limitation of nuclear materials. Additional changes in public demandare also expected due to increasing accessibility of electric vehicles to the average consumer. Thisprojected transition of transportation fuel sources is shown in Figure 1.1.Of course, the growing attention to the impacts of climate change brings evenmore importanceto the effective use of water resources. A steady increase in global surface temperatures sincethe 1970’s has lead to an increased variability of weather patterns, with extreme weather eventsbecoming both more common and more severe (Rahmstorf et al., 2017). California alone hasexperienced 12 drought years since 2000, with the entire US having experienced 5.04 droughtyears on average (N.C.E.I., 2022).These increasing trends of water consumption alongside a more variable climate demonstratehow critical the role of managing water resources as efficiently as possible is. Optimal policiesregarding their use must be the norm, not the exception. Yet this is often not the case, especially indeveloping countries. Thismakes it evermore important to develop transferable, effectivemethodsto determine optimal resource utilization policies. Stochastic dual dynamic programming is onesuch method and is the primary focus of this thesis.3Figure 1.1: Top: Projected Share of New Vehicle Sales. Bottom: Projected Increase inElectricity use in the Transportation Sector (C.E.R., 2018)1.2 Challenges AddressedStochastic dual dynamic programming applications are typically stopped heuristically, e.g. basedon a determined number of iterations and subsequent comparison of lower bounds (stabilization).This means that the definition of a reasonable stopping criterion is left to the user (Fullner &Rebennack, 2021). Recent developments exist which aim to develop deterministic upper boundingtechniques, however these are rather novel and require significant computational effort (Fullner &Rebennack, 2021). The results of this thesiswill assist BCHydro to develop an understanding of theconvergence behaviour of the results from SDDP . This understanding can serve as a foundation tofurther develop deterministic convergence criterion in futuremodel implementations at BCHydro.As well, SDDP problems remain notoriously sensitive to even minimal changes in state-spaceinitial conditions, especially when considering nonlinear functions which impose non-convexitiesonto the future (or expected) value functions. This sensitivity is often referred to as ’algorithmicchaos’ (see Philpott and Guan, 2008 on convergence of SDDP) and can often result in data which4are difficult to interpret (Rougé & Tilmant, 2016). The foundation of SDDP algorithms is the notionthat a future value function can be exactly represented as a piece-wise linear function. This functionis expected to be convex as the nature of a linear representation assumes that a linear ’cut’ wouldbe an appropriate approximation of the function at any given point. These approximations arediscussed in extensive detail in Chapter 2.Figure 1.2: Non-convex Function with Benders Cut (Fullner & Rebennack, 2021)This issue of non-convexity was noted by BC Hydro when attempting to implement a depen-dent random variable into their SDDP model, namely the market price of electricity. Figure 1.2demonstrates how attempting to represent a non-convex function with a linear cut would not aptlyapproximate the function shape. Market prices may be considered a ’random variable’ as theyare indeed stochastic, however they depend heavily on inflow and volume forecasts as well asdemand. Seasonal volume forecasts serve as the primary stochastic variable within the modeldiscussed in this thesis, as the primary focus of this research is to determine if an alternativemodelimplementation of market pricing significantly alters the outcome of the objective function andsubsequently the shape of the piece-wise linear approximation to the future value function. This5new implementation, first suggested by the supervisor of thiswork, Dr. Shawwash, aims to removethe non-convexity issues previously encountered while still considering the markets variance anddependence on demand and inflow forecasts.1.3 Thesis OutlineThis thesis beginswith a thorough and detailed literature review into the development and histori-cal implementation of optimization theory, andmore specifically, stochastic optimization. Chapter2 begins with a brief overview of the history of optimization problems, then followed by a reviewof essential understandings and stochastic dual dynamic programming. Applications in industryare also reviewed from available literature and serve as an example of the efficacy of multi-stageoptimization in the real world.The case study begins in Chapter 3, and provides first an overview of the BC Hydropowersystem and SDDP model implementation. The contents of this section are the primary focus ofthis research. The model presented in §3.1 is the same model in which new approaches regardingmarket price and load demand are being implemented and tested. This model is very dynamicand undergoes changes on an almost daily level, however the subjects discussed within Chapter3 such as treaty / non-treaty storage, water mass balance, etc. remain mostly consistent andwithout change. It is also important to note that BC Hydro’s model was initially developed inAMPL, before being moved to Python and eventually Julia. Julia is a relatively new programminglanguage, as are its packages. This means that, while very quick and efficient, Julia introduces alayer of uncertainty into model performance. Julia is discussed briefly in §3.2.Model results and analysis are provided in Chapter 4. Several assessments are made regardingdifferent performance metrics of the model. Specifically, the existing model (meaning very limitedmarket price dependence on inflow) is compared against a version which attempts to incorporatemarket dependence. The existing model contains market projection information based only on aninflow sequence similar to the year 1972. Ideally, any inflow scenario should be able to referencerelevant market projections specific to that prediction. Some of the metrics considered whencomparing models pertain to: convergence definition, convergence time, computational demands,water value trends, and the shape of Benders cuts.6Chapter 2Background and Literature Review2.1 Theory DevelopmentIn 300BC, optimization played a fundamental role in the development of Euclidean geometry.Today, is resides in almost every facet of modern life. The first consideration of optimizationwas of course by Euclid, the founder of Euclidean Geometry, when he postulated and provedmathematically the shortest distance between two points to be a line (Euclid, 1996). 1900 yearslater in 1615, dynamic programming was invented by Johannes Kepler and applied to a problemwhichwould come to be knownnowas the ’secretaryproblem’ (Ferguson, 1989). However,modernstochastic dual dynamic programmingwas not established until Pereira andPinto publishedMulti-stage stochastic optimization applied to energy planning in 1991. This paper, published in the journalMathematical Programming, described the algorithm for stochastic dual dynamic programmingthat is still largely used today. Theirwork built strongly on that of J.F BendersPartitioning proceduresfor solving mixed variables programming problems, published in 1962. As described by Pereira andPinto, stochastic dual dynamic programming can be interpreted as "Benders cuts in a multistagedecomposition algorithm". As such, Benders cuts and their implementation are heavily discussedin this chapter.The SDDPalgorithmhas beenproven to converge thanks to thework of Linowsky andPhilphottin 2005. On the Convergence of Sampling-Based Decomposition Algorithms for Multistage StochasticProgramswas published in the Journal of Optimization Theory andApplications, wherein Liwoskyand Philpott extend an existing proof by Chen and Powell for the convergence of their ConvergentCutting Plane and Partial Sampling (CUPPS) algorithm to SDDP.72.2 Essential UnderstandingsPrior to developing an algorithm for two-stage or multi-stage stochastic systems, some funda-mentals of linear programming must be addressed. It is assumed that the reader is familiar withsome aspects of optimization theory, however a review of the critical foundations present in thealgorithms discussed in later sections remains important. Since the purpose of this thesis is tobenchmark an SDDP model against itself, the more thorough an understanding of all aspects ofthe model, the more approachable the metrics. This, as well as subsequent sections of this chapter,will not provide proofs; however, many are available pertaining to linear programming. One suchexample is a proof of optimal conditions by Tret’yakov, 2009.The first consideration of fundamental understandings pertains to the concepts of boundednessand feasibility of linear systems. It is understood that the optimal solution always lies at a vertexof its feasible area (Murphy, 2013). While the entire theory of linear programming will not bediscussed here, normal linear maximization problems take the form:maxk∑n=1cnxns.t.k∑n=1ainxn ≤ bi i = 1, . . . , l; n = 1, . . . ,mxn ≥ 0(2.1)where xn, i variables are to be optimized subject to m constraints. Note the objective function is alinear combination of each lth variable. The non-negativity constraint can be circumnavigated bytreating non-positive variables as positive, until they are included in the objective function resultwhere they will return to acting as a negative.2.2.1 FeasibilityThe concept of feasibility is best described as the space (two or higher dimensional) wherein anyvaluewithin that space satisfies all imposed constraints. For a linear twodimensional system, theseconstraints are easily presented graphically. Figure 2.1 demonstrates the feasible region enclosedby three constraints for a linear system of two variables x1 and x2.8Figure 2.1: A Feasible Region Enclosed by Three Constraints (Murphy, 2013)Feasible spaces are considered either semi-infinite (constraints limit either the positive or negativedirection), or finite (fully enclosed by constraints).2.2.2 BoundednessBoundedness is closely related to feasibility, as in order for a linear problem to be consideredbounded, it must not only be feasible, but the objective function must be unable to approach +∞or −∞ for a maximization or minimization objective, respectively. It might be easy to considerthat only finite feasible spaces can be bounded, however, semi-infinite spaces may act as a loweror upper bound depending on the direction of the objective function. Figure 2.2 shows possibleexamples of both bound and unbound problems.2.2.3 OptimalityIn order for an optimal solution to exist for a linear system, the objective must both provide afeasible area as well as be bounded. As mentioned, an optimal solution of linear problems willalways lie at a vertex of constraints. The objective function is also understood to provide an objectivegradient. An objective function will improve fastest along its gradient, and will not improve at allif the feasible space is traversed perpendicular to it. It follows then that a typical pathway foroptimizing a linear program is to first intersect a constraint and then examine how this constraint9Figure 2.2: From Left to Right: a Bounded Finite Feasible Space, a Bounded Semi-Infinite Space, an Unbounded Semi-Infinite Space (Murphy, 2013)aligns with the objective gradient. In fact, moving along the constraints until an optimal solutionis found forms the basis of the Simplex algorithm invented by George Dantzig (Dantzig, 1951).Figure 2.3: Optimality of a Linear System in a Finite Bounded Space (Murphy, 2013)2.2.4 Duality TheorySlightly more complicated than the concepts of feasibility and optimality, is the theory of dualityfor linear optimization problems. In essence, every linear system is considered to have a dual,wherein the dual of a dual problem is the original (primal) problem. This concept is criticallyimportant as it lays the foundation of the stochastic dual dynamic programming algorithm.10From Figure 2.3, it is clear that for any maximization problem, a solution to the objectivefunction which is feasible (but not optimal) provides a lower bound on the optimal value (Murphy,2013). Were the dual problem of that presented in Figure 2.3 considered, any feasible solutionthat is not optimal would instead serve as an upper bound on the optimal value. Consider thegeneral form of a linear maximization problem (2.1). To construct a dual problem from this form,i.e. one which imposes an upper instead of lower bound, the constraints must first be examined,as they are the only components of the system which impose bounds. This section will cover abrief and surface level coverage of duality theory as presented in Murphy, 2013. A more in-depthconsideration of the subject can be found in Vanderbei, 2004.Consider first a weighted sum of the constraints in (2.1), wherein the coefficients of each of thevariables are greater than or equal to those in the objective function. This will impose a bound onthe objective function, so long as the non-negativity condition is true.m∑i=1yik∑n=1ainxj ≤m∑i=1yibiyi ≥ 0 i = 1, . . . , l(2.2)It can then be seen that in order to remain true to the primal problem, the coefficient of eachvariable xn must be at least equal cn. This provides the following set of l constraints:l∑i=1yiain ≥ cin n = 1, . . . ,m (2.3)Since the imposed bound on the objective function is dependent on the selection of scalars yi, it isbest to choose such values which minimize the bounds while still satisfying the constraints. Thus,the dual problem becomes:minl∑i=1yibis.t.l∑i=1yiain ≥ cin i = 1, . . . , l; n = 1, . . . ,myi ≥ 0(2.4)11It is clear from (2.1) and (2.4) that while the primal was a linear system composed of m decisionvariables and l constraints, the dual is now a linear system of l decision variables andm constraints,and has switched from a maximization problem to a minimization. It is critically important toalso understand the concept of normal linear programs, i.e. that the constraints must all be ofthe correct form. This form is ≤ for maximization problems, and ≥ for minimization. Were aconstraint not to be of the correct form, simply multiplying both sides by −1 is enough to bringthe constraint into the normal format. This is known as normalization. Other methods exist to skipnormalization in the construction of a dual problem, however they are not discussed here.The existence of a dual problem can be used to expand upon the definition of optimality brieflydiscussed in §2.2.3. Now, consider the optimal solution to a linear (maximization) system as beingbounded below by a feasible solution to the primal problem, while conversely being boundedabove by another feasible solution to the dual. This is the basis of theWeak Duality Theorem, whichstates that "any feasible solution of the dual is always at least as great as any feasible solution of theprimal" or Z¯ ≤ W¯ (Murphy, 2013). The opposite remains true for a minimization primal. StrongDuality Theorem states that "if the primal problem has any feasible solutions, then the optimal valueof the primal and dual problems is the same" or Z∗ =W∗ (Murphy, 2013). It follows then that thefeasibility and boundedness of the primal is shared by the dual, which carries several implications:• Primal unbounded −→ dual infeasible.• Dual unbounded −→ primal infeasible.• Primal infeasible −→ dual infeasible or unbounded.• Dual infeasible −→ primal infeasible or unbounded.These, as well as knowledge of the formulation of dual problems, provide enough background onthe theory behind linear programming to begin to understand the underlying concepts of stochasticdual dynamic programming, beginning with Benders decomposition.122.3 Decomposable Structures and Benders DecompositionDecomposition is the technique of breaking down an original (non-decomposed) optimizationproblem into n possible decomposed problems, where each decomposed problem is easier to solve.This technique is excellent for simplifying large-scale optimizations such as operational problemsdue to solutions needing to be computationally tractable, as well as exist within a specific timeperiod. Decomposition of these problems means no state variable discretization is necessary, andthus the "curse of dimensionality" encountered in dynamic programming (discussed in §2.4.1) isavoided.However, certain criteria must be met within the optimization problem in order to be decom-posable. Specifically, the problemmust have either complicating variables or complicating constraints.These complicating (sometimes referred to as ’coupling’ or ’linking’) variables are such that, ifconsidered fixed and treated as a parameter, the original problem can be decomposed into a setof subproblems. Similarly, complicating constraints decompose the original problem if they arerelaxed or removed. For example, consider problem (2.5). Notice it can be decomposed due toF6 acting as a complicating constraint. Note that (2.5) as well as subsequent examples within thissection are used with permission from a lecture by Dr. Kazempour of the Technical University ofDenmark (Kazempour, 2022). By removing or relaxing F6, the problem can be decomposed intothree subproblems, one for each set of variables x, y, and z.min A1x1 + A2x2 + A3x3 + B1y1 + B2y2 + C1z1 + C2zz + C3z3s.t. E11x1 + E12x2 + E13x3 ≥ F1E21x1 + E22x2 + E23x3 ≥ F2E31y1 + E32y2 ≥ F3E41z1 + E42z2 + E43z3 ≥ F4E51z1 + E52z2 + E53z3 ≥ F5E61x1 + E62x2 + E63x3 + E64y1 + E65y2 + E66z1 + E67z2 + E68z3 ≥ F6(2.5)Similarly, optimization problem (2.6) contains a constraining variable β. Were β treated as fixed,(2.6) could be decomposed into three subproblems just as (2.5).13min A1x1 + A2x2 + B1y1 + B2y2 + C1z1 + C2zz + D1βs.t. E11x1 + E12x2 + E13β ≥ F1E21x1 + E22x2 + E23β ≥ F2E31y1 + E32y2 + E33β ≥ F3E41z1 + E42z2 + E43β ≥ F4E51z1 + E52z2 + E53β ≥ F5(2.6)The fundamental features of decomposition techniques allow for iterative solutions to both thesingle master problem and the subsequent set of subproblems.Benders decomposition (Benders, 1962) is one such iterative technique to obtain the relationshipbetween the value of the original objective function and the value of the complicating variable.Were the original problem to contain only a complicating constraint, the dual problemwill containa complicating variable and thus the dual can also be solved using Benders decomposition. Asimple example is considered in (2.7):minx,yz = −y− x4s.t. y− x ≤ 5y− 12x ≤ 152y+12x ≤ 352− y+ x ≤ 100 ≤ x ≤ 16y ≥ 0(2.7)By considering x as a complicating variable, it may be removed from the objective problem tocreate subproblem (2.8)14minyiz = −yis.t. yi ≤ 5+ xˆi : πiyi ≤ 152+xˆi2: µiyi ≤ 352− xˆi2: σi−yi ≤ 10− xˆi : γiyi ≥ 0(2.8)where xˆi is the fixed value of the complicating variable x, i is the current Benders iteration, andπ, µ,γ, σ are dual variables obtained through iteration. The notation k will be used to define theset of previous iterations. Therefore if i = 4, then k = 1, 2, 3. The value of the complicatingvariable xmust remain feasible, i.e. it must adhere to the constraint imposed on the original (non-decomposed) problem. Based on the values for dual variables obtained by solving subproblem(2.8), a master problem will obtain the updated value for complicating variable x. In this example,the master problem becomes:minxi ,αi− xi4+ αis.t. 0 ≤ xi ≤ 16αi ≥ ααi ≥ −πk[5+ xi]− µk[152+xi2]− σk[352− xi2]− γk[10− xi]∀k = 1, . . . , i− 1(2.9)where α is referred to as an auxiliary variable and represents the objective function of the subproblem(can also apply to multiple subproblems) within the original objective function. Here, α is a verylarge negative value and is a necessary constraint as at iteration i = 1, k = 0 and thus theminimization would not yet contain a bounding constraint. Additionally, the final constraint is thedual objective function of subproblem (2.8) and is known as a Benders cut. Since the dual variablescome from previous iterations and are parameters, all terms within this function are linear and15thus so are the cuts. Each iteration of the master problem will provide an additional Benders cut.With problem (2.7) now decomposed using Benders decomposition, a simple algorithm can beimplemented to solve the master problem. This algorithm is composed of four steps:1. Set i = 1, xˆ1 = xinitial and lower bound −∞. Subsequent iterations of the master problemalso provide the LB.2. Obtain the values of dual variables π, µ, σ,γ by solving the subproblem. The solution to thesubproblem is the upper bound.3. If |UB - LB| ≤ ϵ, where ϵ is an acceptable level of accuracy, terminate. Otherwise proceedto i = 1+ 1.4. Obtain the updated xi and the corresponding value of αi by solving themaster problem. Thisnew αi is the updated lower bound. Return to Step 1 with xi being the new value for theconstraint variable.It is important to note that the dual objective function of the subproblem represents the lowerbound for the master problem as a consequence of duality theory. Duality theory states that aslong as strong duality holds for a minimization problem, a dual objective function result is lowerthan that of the primal, except at the optimal point where they are equal (Boyd & Vandenberghe,2004). It then follows that the solution to the primal subproblem indeed serves as an upper bound,and vice-versa in the case of a maximization problem.2.4 Optimization using Benders DecompositionTwo-stage problems are simply optimization problems which contain both a first-stage cost, as wellas a future or second-stage cost. Consider an example of a deterministic (non-random) two-stageoptimization problem:16minx1,x2c1x1︸︷︷︸First-stage+ c2x2︸︷︷︸Second-stages.t. A1x1 ≥ b1︸ ︷︷ ︸First-stage constraintE1x1 + A2x2 ≥ b2︸ ︷︷ ︸Complicating constraint(2.10)which can be decomposed to be of the form (assuming a non-negativity constraint):minx1c1x1 + α1(x1)s.t. A1x1 ≥ b1x1 ≥ 0(2.11)minx2c2x2s.t. A2x2 ≥ b2 − E1x1x2 ≥ 0(2.12)where (2.11) is the first-stage problem and (2.12) is the second-stage problem α1(x1). Here, α1(x1)represents the second-stage cost as a function of the first-stage decisions, and is often referred toas the future cost function previously mentioned.Since x1 appears in both the first and second stage, they cannot yet be solved separately. Intheir current state, (2.11) and (2.12) could be solved through discretization of x1 into a set of trialvalues before solving the second-stage problem for each. Intermediate values of the future costfunction α1(x1) can then be obtained via interpolation. Then by using this result, (2.11) can besolved. While feasible for small-scale problems, large-scale optimizations and systemswill becomeexponentially more computationally demanding, and quickly become practicably unsolvable. Forexample, were a problem to contain only 15 decision variables, each with 3 discretized values, atotal of 315 or 14.34 million solutions to the second-stage problemwill need to be completed. Thus,a means to approximate the future cost function α1(x1) by analytical functions instead of discretevalues was developed. This technique is known as dual dynamic programming (DDP).2.4.1 Dual Dynamic ProgrammingConsider once more (2.10) to illustrate the concepts of dual dynamic programming:17min c1x1 + c2x2s.t. A1x1 ≥ b1Ex1 + A2x2 ≥ b2Equation (2.10) can be easily interpreted as a two-stage decision process. The first stage can beconsidered by choosing a feasible trial value for x1, denoted xˆ1. Then, using this value, the optimalsolution for the second stage problem is found. Looking at (2.12) again:min c2x2s.t. A2x2 ≥ b2 − E1 xˆ1Notice the trial value has been treated as a constant and moved to the right-hand side of theconstraints. The objective function remains to minimize the sum of the first-stage and second-stage cost functions. Dynamic programming algorithms (DP) can be used to solve problems suchas (2.10), as well as other stage-based optimization problems. DP defines the first stage as thatshown in (2.11):min c1x1 + α1(x1)s.t. A1x1 ≥ b1Here, the first term represents the immediate cost, and the second the future cost of decision x1. Thiscan be interpreted as the consequences of the first stage decision on the second stage. The futurecost function is defined as (2.12):α1(x1) =min c2x2s.t. A2x2 ≥ b2 − E1x1The future cost function serves to translate the second-stage costs to a function only of the firststage decision x1. These first stage decisions are also known as state variables. As mentioned in 2.4,the DP algorithm attempts to construct the future cost function by discretization and can quickly18result in explosive computation requirements depending on the state space. However, Pereira andPinto, 1991 show that:L(x2,π) = C2x2 + π(b2 − E1x1 − A2x2) + µ(−x2) (2.13)maxπ≥0π(b2 − E1x1)s.t. C2 +−π(A2) ≥ 0(2.14)At the optimal solution, assuming Slater’s condition holds (see Slater, 2013), according to strongduality theorem the future cost function is equal to the maximization of the dual problem (2.14),namely:α1(x1) = π∗(b2 − E1x1) (2.15)where π∗ denotes the value of the dual variable at the optimal point. This shows that it is possibleto formulate a relationship between the future cost function in terms only present in the first-stageproblem. Thus, there exists a linear relationship between x1 and the future cost function so long asthe sensitivity (slope) of π∗ is known. It then follows that in weak duality, (2.15) can be written as:α1(x1) ≥ π(b2 − E1x1) (2.16)where (2.16) is true in any point within the feasible space and not just at the optimal solution. How-ever, it is still not known what the optimal value of π is (π∗). To solve this, assume π1,π2, . . . ,πnis a possible set of values for π wherein an optimal value π∗ exists. It then follows that the futurefunction α1(x1) can be characterized as:19α1(x1) = minααs.t. α ≥ π1(b2 − E1x1)α ≥ π2(b2 − E1x1)...α ≥ πn(b2 − E1x1)(2.17)Each of these constraints can be represented as a cut in the feasible space, allowing for a piece-wise linear approximation of the future cost function without interpolation. This is of course lesscomputationally intense, however, it now seems that not only are the upper and lower bound on πunknown, it is still being discretized. First merging (2.11) and (2.17) since both are minimization,an updated first-stage problem including a piece-wise linear approximation of the future valuefunction is obtained:minα,x1C1x1 + αs.t. A1x1 ≥ b1α ≥ π1(b2 − E1x1)α ≥ π2(b2 − E1x1)...α ≥ πn(b2 − E1x1)(2.18)Expanding on the assumption that the optimal value for the dual variable π∗ exists within thechosen set of π, there is an elegant iterative approach to solve (2.18) known as systematic iterations.The algorithm for systematic iterations is provided in Figure 2.4.Systematic iterations allow for finding the optimal value for π without discretizing π itself.Instead, an initial value for π is selected with an additional constraint being imposed upon it foreach iteration. This system of iterations is the basis of Dual Dynamic Programming and concludesBenders decomposition. A final note: obtaining the optimal solution by Benders decompositionis guaranteed if the objective function of the original (non-decomposed) problem is convex with20Figure 2.4: Algorithm for Generating Possible πi Values for π∗respect to each complicating (linking) variable (Conejo et al., 2011).2.4.2 Stochastic Dual Dynamic ProgrammingStochastic (random) two-stage problems are very similar to the example provided in §2.4.1. Thedifference is found in the possible outcomes of the second-stage results; in deterministic two-stageoptimization, there is only one possible second stage scenario. That is to say, one outcome with100% certainty. Second-stage stochastic decisions have a distribution of outcomes each with adistinct likelihood or probability of occurring. There is no practical limitation on the possiblenumber of outcomes, ranging from 1, 2, . . . , n. Stochastic second-stage scenarios are generatedbased on available forecasts in the first-stage results (or day-ahead stage). The general form of atwo-stage stochastic optimization problem is:minxω ,yω∑ωφω [ f (xω) + g(yω)]s.t. h(xω) ≤ α ∀ωk(xω) + l(yω) ≤ bω ∀ωxω = xω˜ ∀ω, ∀ω˜(2.19)21where φω corresponds to the probability of scenario ω. The set of first-stage variables arescenario independent, meaning they do not depend on the probabilities of each second-stageoutcome. This is ensured by imposing the non-anticipatory condition xω = xω˜. This constraint isstating that for all sets of scenariosω and ω˜, the optimal decision for first-stage variable xwill be thesame. It is referred to as the non-anticipatory condition as it is imposed to prevent the anticipation offuture information. However, the second-stage variables are of course scenario dependent. Sincethere is no certain minimum value for a stochastic problem, (2.19) is minimizing the expected value,a concept frequently utilized from probability theory.The example presented in (2.19) is perfect for a simple two-stage problem. However, oftenreal-world applications of this technique are multi-stage, requiring several optimizations overnumerous time periods. Before generalizing into a multi-stage framework, consider first a three-stage stochastic problem:Figure 2.5: Three-stage Stochastic Scenario TreeFor the hypothetical problem presented in Figure 2.5, future value functions still need tobe solved, however there now exist dependencies where third-stage decisions depend on thedecisions of the second-stage. Utilizing the general form of (2.19), a three-stage stochastic problemis presented as:22minxω ,yω∑ωφω[f (xS1ω ) + g(xS2ω ) + h(xS3ω )]s.t. i(xS1ω ) ≤ α ∀ωl(xS1ω ) +m(xS2ω ) + n(xS3ω ) ≤ bω ∀ωxS1ω = xS1ω˜ ∀ω, ∀ω˜xS2ω = xS2ω˜ if ΘS2ω = ΘS2ω˜ ∀ω, ∀ω˜(2.20)where the final constraint is imposing a non-anticipatory condition on the second-stage deci-sions in addition to those imposed on the first-stage in (2.19). Here, Θω represents the randomvariables related to the corresponding stage (S2). Consider two hypothetical scenarios ω1 and ω2:ω1 =ΘS21 = 10ΘS31 = 20ω2 =ΘS22 = 10ΘS32 = 25Since in S2 the random variables ΘS21 and ΘS22 are equivalent, the corresponding S2 decisionsxS21 and xS22 are also equivalent.2.4.3 Multi-stage Stochastic Problems and Nested Benders DecompositionThe two and three-stage stochastic problems presented in §2.4.2 can be further generalized intomulti-stage stochastic problems of the form:minxt,ω∑ωφω∑tf (xt,ω)s.t. g(xt,ω) ≤ αtω ∀t, ∀ωxt,ω = xt,ω˜ i f Θt˜,ω = Θt˜,ω˜ ∀t˜ ≤ t, ∀t, ∀ω, ∀ω˜(2.21)where the new variable t represents the index of stages and t˜ is the index of previous andcurrent stages. Thus, if t = 5, then t˜ is the set of all stages t1 . . . t4. The generalized form presentedin (2.21) can be solved by utilizing a technique known as Nested Benders Decomposition.Nested Benders Decomposition solves multi-stage stochastic problems by breaking down thesystem into a nested series of two-stage problems. Figure 2.6 shows a four-stage stochastic system23Figure 2.6: A Typical Stochastic Decision Tree (Murphy, 2013)which can be represented by 13 two-stage problems, 12 of which are nested (Figure 2.7). Circledin red in Figure 2.7 is the master problem for all subproblems. As was the case with non-nestedBenders decomposition, an iterative technique will be used to solve for the optimal solution.However, each iteration will now need to solve several two-stage problems.Figure 2.7: Nested Two-stage Problems (Kazempour, 2022)Consider an example from Murphy, 2013, wherein nested Benders decomposition is used tosolve a three-stage stochastic problem, presented in Figure 2.8. The second and third stage ofthis problem contain three (B, C, D) and six (E, F, G, H, I, J) nodes, respectively. As previously24Figure 2.8: A Three-stage Stochastic Problem (Murphy, 2013)mentioned, each iteration will now have to solve several nodal subproblems. The first iterationwill be in the forward direction, meaning the problem will be solved starting from the first stageas follows:1. Solve the master problem A. Stage two nodes are now nodal subproblems. The notation A1refers to node A, iteration 1.minxA1,αA1f (xA1) + αA1s.t. g(xA1) ≤ bαA1 ≥ αmin(2.22)αA1 represents the future cost from nodes B, C, and D.2. Solve nodes B, C, and D. These nodes are now master problems for nodes E to J, whileremaining subproblems to node A. Consider node B:minxB1,αB1,xA1i(xB1) + αB1s.t. j(xB1) ≤ c− k(xA1)αA1 ≥ αminxA1 = xA1,fixed : γA1→B1(2.23)25αB2 represents the future cost from nodes E, F. Note the addition of local constraints as afunction of the master problem decision variable xA1. γ is the corresponding dual variablefrom decision xA1.3. Solve nodes E to J. Consider E, a subproblem of B:minxE1,xB1l(xE1)s.t. m(xE1) ≤ d− n(xB1)xB1 = xB1,fixed : γB1→E1(2.24)There is no future cost function (α) here as it is the final node associated with the problem.A critical assumption is made here, specifically that node E (or F, G, etc.) does not depend on thedecisions of node A. Thus it is assumed that subproblems only depend on the decision variables oftheir associated master problem. The second iteration will now move in the backwards direction,solving from nodes E to J back to node A. This iteration will introduce Benders cuts.1. Since E to J have already been solved from the previous iteration, B, C, and D will beconsidered next, utilizing associated dual variables from subproblem nodes E to J . Considernode B once more:minxB2,αB2,xA1i(xB2) + αB2s.t. j(xB2) ≤ c− k(xA1)αB2 ≥ αminxA1 = xA1,fixed : γA1→B2αB2 ≥[l(xE1) + o(xF1)]+[γB1→E1 + γB1→F1] (xB2 − xB1)(2.25)The final constraint in (2.25) is the Benders cut that will be used to aid in the determinationof a piece-wise linear approximation of the future cost function of the master problem.2. Finally, solve the master problem A, including Benders cuts.26minxA2,αA2f (xA2) + αA2s.t. g(xA2) ≤ bαA2 ≥ αminαA2 ≥[i(xB2) + αB2 + p(xC2) + αC2 + q(xD2) + αD2]+[γA1→B2 + γA1→C2 + γA1→D2] (xA2 − xA1)(2.26)Several iterations, both forward and backward, will need to be completed until the future costfunctions converge. However, not all nodes necessarily need to reach a convergence criterion. Sincethe Benders cut in (2.26) contains the future cost function of all subnodal problems, convergenceof the root node is usually sufficient (Kazempour, 2022). However, depending on context of themodel, it may be pertinent to check convergence at each node.2.4.4 Convergence Proof and Criterion of SDDPModelsProblems such as the example presented in §2.4.3 are known as ’multi-stage stochastic linearprograms with recourse’, and have had vigorous mathematical proofs conducted detailing theirconvergence (Philpott & Guan, 2008). One of the earliest statements regarding SDDP convergenceis found in a PhD dissertation by Donohue, 1996, in which it is stated, without a formal proof,that convergence of this algorithm follows from the finite convergence of the nested Bendersdecompositionalgorithm, as each scenarios optimality cuts are re-sampledat each iteration. Formalmathematical proofs began with the Cutting-Plan and Partial Sampling algorithm (CUPPS) proofby Chen and Powell, 1999, which was then expanded upon by Linowsky and Philpott in 2005 toencompass the SDDP algorithm (Linowsky & Philpott, 2005).However, the proofs provided in both Chen and Powell, 1999 and Linowsky and Philpott,2005 did not provide an explicit assumption regarding the independence of random variablesand convergent subsequences of each iteration. Philpott and Guan argue that the absence of thisassumption "weakens the analysis in these papers, and leaves open the question of convergencein general" (Philpott & Guan, 2008). Their 2008 proof provides what is described as a "simplerproof of almost sure convergence" for several sampling algorithms, SDDP included, which does27not require the assumption previously mentioned regarding independence of random variables.Instead, the proof follows the arguments presented in Donohue, 1996, now explicitly stated asassumptions. The derivation of their proof lies beyond the scope of this thesis, however, the workof Philpott and Guan serves as the foundation for confidently setting convergence criterion for anySDDP model contained herein.Although SDDP models have been proven to converge under pertinent assumptions, the cri-terion for determining convergence may vary from model to model, depending on context. Theconvergence criterion provided in §2.4.3 wherein future value functions are compared betweennodes is one such option. An alternative option discussed in §3.1 relies on Monte Carlo Sampling(MCS), which to the knowledge of Guan et al., 2018, had not been implemented into an SDDPmodel prior to their program development.2.5 Applications from LiteratureStochastic optimization using dual dynamic programming is an excellent approach to solve large-scale multireservoir hydroelectric problems. While the objective of each model may vary based onthe needs of the system, such as unit commitment, dispatch, scheduling, or energy exports, SDDPhas frequently been utilized to provide computationally tractable solutions to each.2.5.1 Optimal Scheduling of a Hydrothermal Generating System by Pereira & PintoPereira and Pinto, 1991 proposed an SDDP approach to solving the optimal scheduling of a hy-drothermal generating systemwhich, for each stage of a specified planning horizon and associatedstate variables, could provide generation targets for each generating facility. Their model aimedto minimize the total expected cost of operation, specifically fuel costs and penalties for failure toprovide adequate load supply. Pereira and Pinto note a very common occurrence in these models,specifically the complexity imposed by state variables at every stage and the future consequencesthey impact by linking previous decisions to those in the future through, for example, reservoirstorage at each stage. As they put it, the model must balance between management of resourcesas well as revenue.28Aswith allmodels of this nature, the stochastic variables arise from imperfect forecasts of inflowsequences and/or market costs. In Pereira and Pinto’s case, the nonlinear generation costs and theexpressed form of outflow and head in hydroelectric production result in an nonlinear objectivefunction which is also nonseparable due to the worth of energy generated in a specific plant notbeing a direct function of the plants state alone. This goes to show that, while the theory examplesprovided in this thesis apply to both separable (or decomposable) and linear objectives, SDDP isstill very much applicable to nonlinear inseparable problems, albeit inherently more complex toimplement. The requirement of problem convexity remains true for the nonlinear case.αt(xt) = Eat|xt{minutct(ut) + βαt+1(xt+1)}s.t. xt+1 = ft(xt, at, ut)gt+1(xt+1) ≥ 0ht(ut) ≥ 0, ∀t = 1, . . . , T − 1, T ∀xt(2.27)wheret: indexes the stages for T planning horizon,xt: state vector at the beginning of state t,αt(xt): expected value of operation cost from state xt,at|xt: probability distribution of inflow vector at conditioned by state xt,E{·}: expected value,ut: decision vector for the stage,ct(·): immediate cost associated with decision ut,β: discount factor,xt+1: state transition equation,gt+1(xt+1): constraints on the state vector,29ht(ut): constraints on the decision vector.The system is further explained in detail in Pereira and Pinto, 1991.The objective function for this model assumes the inflow volumes are known at the beginningof each stage. This assumption is shared by the system discussed in §3.1 and plays a significant rolein model execution. The objective is applied to a system composed of 39 hydroelectric plants (22 ofwhich contained reservoirs, and 17 being run-of-the-river). Initial state variables were always setto 50% of the storage capacities. Thermal generation costs and energy demands were also fixed.Inflows were represented by two independent random vectors with two possible realizationsat each stage with the first stage is assumed to be known. This provided 512 possible inflowsequences. In total, the model considered 1023 subproblems with 104,346 variables and 40,920constraints. The problemwas solved utilizing a network flow algorithm developed by Kenningtonand Helgason, 1984 which took advantage of the similarities between the branching SDDP systemstructure and network flows.Pereira and Pinto reached an optimal solution in ’five iterations’ (10 iterations if the definitionof an iteration provided in §2.4.3 is considered instead). The lower-bound was obtained from thesolution to (2.27), whereas the upper-bound was estimated from Monte Carlo simulation (MCS)results for all stages and scenarios:z¯ = c1x1 +1nn∑i=1zi (2.28)where n is sample size and zi is the total cost for one Monte Carlo run, specificallyzi =T∑t=1ct xˆti. (2.29)Pereira and Pinto utilized a confidence interval as a convergence criterion through the calculationof standard deviation viaσz =√1n2n∑i=1(z¯− zi)2. (2.30)30Figure 2.9: Expected Operations Cost (Pereira & Pinto, 1991)Computation time to solve this problem was 37 minutes utilizing 1990’s technology with a CPUclock speed of 20 MHz. The limited inflow probabilities (only two per stage) coupled with loadand cost values being fixed, greatly simplified this large-scale problem down to a tractable leveleven in the 1990’s.2.5.2 Optimal Operation in Deregulated Electricity Markets by Iliadis et al.Similar to themodel developed byPereira andPinto, 1991, Iliadis et al., 2008 presented an optimiza-tion problemwhich considered not only stochastic inflows, but energy prices as well. Randomnesswas considered for both spot and forward prices with the objective of maximizing revenue bysimultaneously determining generation and selling / purchasing decisions for current and futuremarkets. This more closely represents the system discussed in §3.1, however Iliadis et al. considera market dominated by monopolies and oligopolies which shifted the interests of companies fromcost minimization (as was seen in Pereira and Pinto, 1991) to revenue maximization.Inflows for this system were modelled as a variable lag auto-regressive AR(p) model. Statevariables pertained to reservoir storage levels and inflows from previous months. Since theobjective of this model is to maximize revenue, the sum of future benefit (future value function) aswell as the immediate benefit must be considered. Iliadis et al. refer to these functions as future31benefit (FBF) and immediate benefit (IBF). It is assumed that this power company is a ’price taker’,meaning it cannot influence the market prices by selling.Figure 2.10: Optimal Point of FBF and IBF Sum (Iliadis et al., 2008)Prices were considered to be auto-correlated, meaning the current stage price carried infor-mation about the possible outcomes of future prices, and as such were treated as state variables.This imposed significant issues on the FBF as it removed concavity for some state variables, andthus removed the ability to predict the FBF through linear Benders cuts. To solve this, Iliadis et al.removed cut sharing across price states, meaning the FBF was constructed for each price state ateach stage. The objective function thus became:βt(vt, qt−1,πt) = Eqt,πt+1|qt−1,πt max{πtet −∑jc(j)gt(j) + βt+1(vt+1, qt,πt+1)}(2.31)whereqt: inflow to the reservoir during period t,πt: price for stage t,vt: volume in storage at the beginning of time t,32βt+1: future benefits function,E{·}: expected value.Constraints were not provided in Iliadis et al., 2008, however they may be possible to deducethrough a thorough review of the paper.This model was applied to a hydropower system in France consisting of 35 hydropower plantswith 15 reservoirs. Inflows and prices were randomly combined from 50 possible scenarios withfive weekly price blocks. Green market prices were assumed to be a guaranteed known price. Aweekly time stepwas usedwith spot price scenarios aggregated intoweekly price blocks to captureintra-week decisions. This was done to level out the volatility of pricing. Additionally, these blockswere aggregated in such a way that hours which typically have similar behaviour were groupedtogether. The means of developing these price blocks was not disclosed and simply referred toan ’in-house’ model. This, along with not providing model constraints, was likely due to thedevelopment of the model being secretive as to not provide benefit to the competing hydropowercompanies. Convergence criteria was also not discussed. Figure 2.11 provides the model resultsrun over a 12 month planning horizon.Figure 2.11: Optimal (12Month) Hydropower Generation Results (Iliadis et al., 2008)33Chapter 3Case Study3.1 Model of the BC Hydropower SystemBritish Columbia, Canada (BC) is home to the Columbia River, the largest river in the PacificNorthwest of North America. At over 2,000km long, the Columbia River is facilitates over 60dams spanning across both Canada and the United States. In 2012, the Columbia River accountedfor one third of all hydroelectric potential of the United States and over 40% of all hydroelectricpower production (Lillis, 2014). It follows that the importance of optimal operations of Canadianhydroelectric projects cannot be understated. In 2017 Dr. Shawwash, alongside Dr. Guan andDr. Abdalla, published a paper detailing the development of two novel SDDP models of BC’smultireservoir hydropower system (Guan et al., 2018). The models developed in their paper focuson generating water value functions for operations planning at BC Hydro.BCHydro is a provincial electric utility andCrown corporation serving 95%of BC’s population.Primary hydropower generation for BCHydro comes fromfive facilities, namely the GoronMerrittShrum (GMS, 2820MW), Peace Canyon (PCN, 700MW), Mica (MCA, 2806MW), Revelstoke (REV,2518MW), and Hugh Keenleyside (ARD, 185MW) dams. The location of these facilities, as wellas the locations of the Peace River and Upper Columbia River, are presented in Figure 3.1. Guanet al., 2018 considers these facilities in operations planning to make decisions regarding reservoirreleases to meet forecasted domestic electricity demands as well as maximize total benefits fromenergy trade in US and Alberta markets. Of specific importance to this thesis is the considerationof the value of water stored in reservoirs over a multiyear planning horizon. Reservoir releasepolicies must consider stochastic inflows, prices, and loads, as well as complywith Columbia River34Figure 3.1: Network of Primary BCHydro Hydroelectric Facilities (Guan et al., 2018)Treaty (CRT) requirements and all related agreements (Guan et al., 2018).3.1.1 Model DevelopmentThe nature of a multireservoir system such as that operated by BCHydro necessitates the develop-ment of a model which can incorporate as many aspects of the system as possible. In this case, twoadditional models feed into the SDDP optimization, namely an Inflow Model and a CRT Model.As with any hydropower facility, operations depend heavily on the inflow forecasts for that par-ticular facility. Additionally, the special nature of operations under the CRT require specific floodcontrol measures and storage requirements to exist at specific facilities. Random variables (RV) areintroduced by two stochastic models consisting of seasonal volume forecasts and monthly inflows(together known as the InflowModel), whereas limitations imposed by the CRT are known. Thesetwo models are briefly discussed in this section. For more thorough breakdowns of model com-ponents, see the original paper by Guan et al., 2018. The networking of all implemented models isshown in Figure 3.2.35Figure 3.2: Relationships Between Implemented Models (Guan et al., 2018)The objective of this model is to determine optimal operation of the five main hydroelectricfacilities at monthly stages, dependent on the storage and hydrological states of the system andrealization of the stochastic models. This can mathematically be represented by the system:Ct(Xk,m,t,ωt)D.V. = max ∑k∈K,m∈MBk,m,t(Xk,m,t) + Bdsbt (Ut)−Ct(P.V.) + γEωt+1∈Ωt+1 [Ct+1 (Xk,m,t+1,ωt+1)](3.1)whereCt(Xk,m,t,ωt) is the objective function at stage t,D.V. is the list of decision variables,P.V. is the list of penalty variables,Bk,m,t is the concave piecewise linear function to calculate benefit from net export Xk,m,t,Bdsbt is a function that calculates downstream benefits as a function of release Ut,36γ is the fixed discounting factor for discounting of water value (roughly 0.995),ωt is the vector of random variables sampled from Ωt where ω1 is known and fixed,Ct is the linear function to calculate the cost from penalty variables,Eωt+1∈Ωt+1 [Ct+1 (Xk,m,t+1,ωt+1)] is the expected value of Ct+1(Xk,m,t+1,ωt+1).Four sets of constraints are imposed on (3.1), each derived from modelling specific storagerequirements, namely flood control, treaty and non-treaty operations, and water/mass balance.Inflowcalculationsmust also be discussed to understand the selection of seasonal volume forecasts.InflowModelThe development of the inflow model utilizes 60 years of historical monthly inflow data as well asseasonal volume forecasts. This data is processed using five statistical methods: linear regression,an autoregressivemodel of order 1 (AR1), shifted log-transformation, principal component analysis(PCA), and probability distribution approximation (Guan et al., 2018). Seasonal volume forecastsbeginning at each month are generated alongside the monthly inflows for Peace River and upperColumbia River prior to being distributed into corresponding river sub-catchments proportionalto historical allocations. The flow of data processing in the Inflow Model is presented in Figure3.3.Peace River and upper Columbia River inflows typically pick up during the freshet period. Thisperiod begins in April and peaks in either June or July, before significantly decreasing followingmid-August. Inflows during this time primarily source from snow melt and run-off from precipi-tation. Run-off is typically higher during the warmer summermonths surrounding August. Theseinflow volumes are highly dependent on the quantity of snow pack that accumulated within thecatchment during the colder months of November to April. Information pertaining to the precisesize of the snow pack is available to BC Hydro usually by late December or early January. AfterApril, this size confidently knownwith a high accuracy. The subsequent seasonal volume forecastsfor both the Peace and upper Columbia River are used for decision making purposes in operationof reservoirs beginning in early January. Subsequent inflows that do not come from snowmelt aremainly from run-off due to precipitation. Consequently, inflows from August to November tend37to be large and have significant variations. These variations taper off come December to March, asmost of the precipitation has begun to accumulate as snow pack.Figure 3.3: Use of Statistical Methods in Inflow RV Development (Guan et al., 2018)From January to April, seasonal volume forecasts are generated using linear regression. Theresiduals of this regression are modelled by the AR(1) process due to a strong ’lag-1’ correlation.Recall that a lag-1 correlation refers to a relationship between two points that are one time periodapart. White noise which may arise from the AR(1) model for March and April is typicallyignored due to the strong residual correlations found in February to April. As such, the shiftedlog-transformation is not required due to the noise being normally distributed. PCA is thenused to calculate the principal components from the white noise, which are used to calculatespatially correlated white noise from two independent normally distributed random variables ineach month. Lastly, a discrete probability distribution is generated with a number of distinctvalues and probabilities. These are used to approximate the normal distribution and sample thetwo independent random variables in each month.From this,monthly inflowsaregeneratedusing linear regression. These inflowsare conditionedon the seasonal forecasts. To properly represent the stochastic behaviour of inflows that originatefrom snow melt, seasonal volume forecasts are used as hydrological state variables. With threesampled distinct values, in total there are nine outcomes of forecasts (nine for both the Peace Riverand upper Columbia River) in January, and 81 from February to April.However, for April to November, the AR(1) process models the residuals owing to their stronglag-1 correlation. The generated white noise is not ignored here and is instead transformed by theshifted log-transformation into a normal distribution. This transformed white noise can then beassessed by PCA. Since inflows from April to August source mainly from snow melt, and those38from September to November from run-off due to precipitation, the AR(1), log-transformation,and PCA for these inflow periods must be developed separately. The probability distributionapproximation is conducted as previously mentioned for generating seasonal volume forecasts.Thus the inflowmodel also generates nine outcomes for both the Peace and upper Columbia Riversfrom April to November for each outcome. From April to November, the model generates 729 to3.5 billion scenario outcomes of inflows.Flood ControlDuring the freshet period fromOctober to April, ensuring adequate storage availability to mitigateflood risk is crucial. Specifically, these storage constraints are imposed on MCA and ARD. Theseconstraints are modelled as:Aj,t = aj,t, j ∈ MCA, ARD ∀t ∈ October to DecemberAj,t = max{min[Aj,t(Ftdt ), Aˆj,t], Aˇj,t}∀t ∈ January to AprilVj,t ≤ Vˆj − Aj,t ∀t(3.2)whereAj,t is the storage space for flood control at stage t,aj,t is fixed storage space,Aj,t is a linear regression equation developed using Storage Reservation Diagrams (SRD),Ftdt is the seasonal volume forecast,Aˆj,t and Aˇj,t are the maximum and minimum storage limits respectively,Vj,t is the storage of reservoir j at stage t,Vˆj is reservoir storage capacity.Treaty StorageGoverned by the Columbia River Treaty, the total treaty storage accounts in MCA and ARD at theend of each month must meet the total treaty storage that is set for the beginning of the month.39Thus the total treaty storage depends on that of the previous month, inflow volume forecasts, aswell as monthly inflow volume forecasts and inflow from both the upper Columbia River and TheDalles. A breakdown of the modeled reservoirs and their respective storage accounts is shown inFigure 3.4.Figure 3.4: Modeled Reservoirs and Storage Accounts (Guan et al., 2018)These storage requirements serve to ensure the US is protected from any negative downstreamconsequences of unmitigated release policies. The treaty constraints can be represented as:Ytott = min{max[Yt(Ytott−1,Qt, Ft), Yˆtott ] , Yˇtott } ∀t ∈ January to JuneYtott = min{max[Yt(Ytott−1,Qt), Yˆtott ] , Yˇtott } ∀t ∈ July to DecemberUARDt =1dtYtott−1 +Qt −1dtYtott , ∀t(3.3)whereYtott is the total treaty storage at stage t,Yt is a regression equation, Qt corresponds to monthly inflows,Ft corresponds to forecasts of inflow volumes,Yˆtott and Ytott are the maximum and minimum limits for total treaty storage,UARDt is the ARD treaty release,dt corresponds to the number of days in t.40However, in real operation, the treaty storage inMCA and ARD is predetermined. Flex storageaccounts are used to shift the amount of water stored between MCA and ARD treaty amounts.Thus the constraints can be more adequately represented as:YMCAt = YMCAt−1 + dtQMCAt −1hday ∑k∈Khk,tUMCAk,t ∀tYMCAt +YARDt = Ytot ∀tYˇMCAt ≤ YMCAt ≤ YˆMCAt ∀tYˇARDt ≤ YARDt ≤ YˆARDt ∀t(3.4)wherehday is the number of hours in a day,hk,t is the number of hours in price-load blocks k.Non-Treaty StorageNon-treaty storage accounts for MCA follow conditions outlined by the Non-treaty Storage Agree-ment between BC Hydro and the Bonneville Power Administration. The US and BC non-treaty(BCN) accounts have identical storage capacities and it is assumed in operations planning thatUS non-treaty (USN) storage and release will match BC’s non-treaty storage and release. Theconstraints imposed by operations of non-treaty accounts can be represented by:YBCNt = YUSNt ∀tUBCNt = UUSNt ∀tYBCNt = YBCNt−1 − dtUBCNt ∀tYˇBCNt ≤ YBCNt ≤ YˆBCNt ∀tUBCNt +UBCN+t ≥ UˇBCNt ∀tUBCNt −UBCN-t ≤ UˆBCNt ∀t(3.5)whereUBCNt ,UUSNt are BC and US non-treaty releases,41UBCN+t is a non-negative penalty variable for not meeting release requirements,UBCN-t is a non-negative penalty variable for exceeding release limit.Water Mass BalanceMass balance is a fundamental building block of hydrology theory andof course plays an importantrole in defining some of the constraints which need be imposed on this model. The water massbalance equation for all storage facilities in this system is represented as:V t−1 +1hday ∑k∈Khk,t (Qt −GRk,t −GSt) = V t (3.6)whereV t−1 is the vector of initial storage (states) for the five reservoirs,Qt is the vector of monthly inflows,G is the matrix which defines the hydraulic flow path,Rk,t is the vector of turbine releases in price-load block k,St is the vector of spillage amounts,V t is the vector of storage for all five reservoirs.For GMS, MCA, and ARD, the storage is constrained by maximum and minimum limits,whereas PCN and REV are fixed at essentially full storage due to their small capacities. Turbinereleases and spillage are constrained by maximum limits as functions of averages of initial storagestates and decision variables for GMS, MCA, and ARD. However, these are fixed for PCN andREV. Penalty variables may be implemented to relax these storage requirements. Turbine releasestake into account a fixed monthly maintenance schedule. Plants adhere to minimum releaserequirements which consider the total of turbine release and spillage. Limits on release for plantsGMS and PCN are influenced by ice flow limits of the Peace River downstream of Site C in thewinter months.42Energy transmission must also adhere to a load and resource balance as well as transmissionlimits for energy export. These limitations are presented as:∑j∈JPj,k,t + Pshk,t + Pfixt − ∑m∈MXk,m,t − Xfixt + L+k,t ≥ Lk,t ∀k ∈ KXˇk,m,t ≤ Xk,m,t + Xfixt ≤ Xˆk,m,t∀k ∈ K(3.7)wherePk,m,t is the power generation of plant j in price-load block k,Pshk,t is the power generation of small hydro plants owned and operated by BC HYdro,Pfixt is monthly energy from other energy resources such as IPP’s and is assumed to be fixed,Xfixt is pre-scheduled energy export and is assumed to be fixed,L+k,t is a penalty variable for not meeting the load requirements,Lk,t is the system load,Xˇk,m,t and Xˆk,m,t are transmission limits for energy export.3.2 Model Implementation in JuliaThe development of an SDDP model which adequately represents the interactions and utilizationof BC Hydro’s reservoir systems began in AMPL before migrating to Python and finally Julia.Original AMPL implementations saved the linear programming (LP) model at each stage beforeutilizing a solver (e.g. CPLEX) to recall and solve the LP model at each iteration. This of courseproved computationally heavy and inefficient for large-scale models.This issue was overcome by transferring the model to a Python environment. The nature ofPython, specifically its high-level syntax and built-in data structures allowed the step of savingthe LP model to be removed. While the Python model performed better, interpretive languagessuch as Python rely on instructions being executed by an external program rather than the targetmachine. This means that Python can provide an overall performance increase from AMPL, butfalls behindwhen compared to compiler-based languages such asC++and Julia. This thesis utilizes43an application of SDDP algorithms in the Julia language to achieve further reduced computationaltime. It is important to note that the Julia model was fully functional prior to the developmentof this thesis, and the work contained herein serves only to determine if specific changes to theexisting model could achieve better performance based on specific metrics.3.2.1 SDDP AlgorithmRecall from §2.4.2 that SDDP is an efficient approach to solvingmultistage stochastic programmingproblems of the formZt(xt−1,ωt) = maxxt(c⊤t xt +Eωt+1∈Ωt+1 [Zt+1(xt,ωt+1)])s.t. Atxt = ωt + Bt−1xt−1(3.8)whereZt(xt−1,ωt) is the optimal objective function value,xt is the vector of state variables,c⊤t is the vector of coefficients of xt in the objective function,At is the matrix of coefficients of xt in the constraints,Bt−1 is the vector of coefficients of xt−1,Eωt+1∈Ωt+1 [Zt+1(xt,ωt+1)] is the expected value of Zt+1(xt,ωt+1), andωt is the vector of random variables sampled from Ωt.Note that the notation used to represent an SDDP system has changed from previous sectionsand now will be represented in terms of (3.8) to be in line with the notation used at BC Hydro.For feasibility, it is assumed that the expected value of the objective function at t+ 1 is a concavefunction of the state variables such that linear Benders cuts are constructed as a function of thestate variables. The minimum of these Benders cuts represents the outer approximation of thewater value function as a piecewise linear function, shown in Figure 3.5.44Figure 3.5: Benders Cuts and water value ApproximationTo solve (3.8) at stage t, the SDDP model isZt(xt−1,ωt) = maxxt(c⊤t xt + θt)s.t. Atxt = ωt + Bt−1xt−1θt ≤ αt,n + β⊤t,nxt, n ∈ N, t < Tθt = Eωt+1∈Ωt+1 [Zt+1(xt,ωt+1)] t = T(3.9)whereθt is the water value,αt,n + β⊤t,nxt are Benders cuts where αt.n is the intercept of Benders cut n at stage t, and βt,n is thevector of slopes of xt of Benders cut n at stage t.Benders cuts are generated in passes, wherein each pass consists of a forward and backwards pass.The algorithms for both the forward and backwards passes are provided in pseudo-code below.Algorithm 1 (forward pass) provides optimal trajectories of the state variables (the vector xt)given the initial state x0. The maxima and minima of each vector for each stage t of the optimaltrajectories form the optimal ’corridor’ for xt.45for t = 1, . . . , T − 1 dofor s = 1, . . . , S dosample ωt from Ωt where ωt is known and fixed.solve stage t with xt−1(s) and ωt where x0(s) is known and fixedrecord xt(s)endendAlgorithm 1: Forward Passfor t = T, . . . , 2 dofor s = 1, . . . , S dofor ωt ∈ Ωt dosolve stage t with xt−1(s)record Zt(xt−1,ωt) and λt(ωt) as the Lagrangian multipliers of the equalityconstraints (not including the final water value function) in stage tendcalculate αt−1,n+s and βt−1,n+s using Zt(xt−1,ωt) and λt(ωt) for ωt ∈ Ωt asαt−1,n+s = Eωt∈Ωt[Zt(xt−1(s),ωt)− λ⊤t (ωt)Bt−1xt−1(s)]andβt−1,n+s = Eωt∈Ωt[B⊤t−1λt(ωt)]and add the Benders cut to the stage t− 1 problemendendAlgorithm 2: Backward Pass3.2.2 Model ConvergenceGiven sufficient iterations, the Benders cuts provide a reliable approximation to the true watervalue functions, and need only be tested for convergence. Techniques developed by Philpott andGuan, 2008 provide convergence criterion that is directly applicable to SDDP algorithms of thetype used in the development of this water value model. These criterion build off the formalproof of convergence for CUPPS provided by Chen and Powell, 1999 which was then extended tocover other multi-stage sampling algorithms such as SDDP by Linowsky and Philpott, 2005. Thesecriterion require:1. Random quantities appear only on the right-hand side of the linear constraints in each staget,2. The set Ωt of random outcomes in each stage t = 2, 3, . . . , T is discrete and finite (Ωt ={ωti|i = 1, . . . , qt < ∞},3. Random quantities in different stages are independent,464. The feasible region of the linear program in each stage is non-empty and bounded.These criterion hold true for the Julia formulation of the SDDP algorithm and thus the proofprovided by Philpott and Guan, 2008 holds for this case. Thus, the results of the SDDP systemmust only be tested to see if convergence has yet been obtained provided a certain number ofiterations. The converge criterion for each node is assessed in accordance to the steps:1. Run algorithm 1 and 2 to generate Benders cuts and solve the t = 1 problem. RecordZ1(x0,ω1),2. Runadditional passeswith thegeneratedBenders cuts to calculate the total return∑Tt=1 c⊤t xt+θt for each scenario. For each scenario, calculate the average total return z which is an esti-mate of the lower bound of Z1(x0,ω1). Note that this step is the simulation process used toassess the WV functions,3. Calculate a 95% confidence interval of z as (z− 1.96σ, z+ 1.96σ),4. If Z1(x0,ω1) is within this interval, go to step 5, otherwise return to step 1,5. Check convergence criteria. If satisfied, terminate. Otherwise, runnewpasses to addBenderscuts.A similar test is used to determine overall convergence whereby the results of the simulation areused to determine overall water value (θ) at each stage. Iterations are incremented in stages of500 beginning from 500 and ending at 3,000. Subsequent iterations are increased by 1,000. Theresults of the final stage of simulation (27 years for the market data incorporated model, 49 yearsfor the unalteredmodel) are compared to the results of the previous iteration run. If the results aresufficiently close at all stages, say within some percentage of error ε, then the model is consideredto have converged. This criterion is simple, yet comes with some downsides. Often times, nodes(stages) will approach an optima from different directions, depending on the number of iterations.This is sometimes seen when using Gradient Descent Optimization and can be managed througha technique known as aggregated momentum (Lucas et al., 2019). Interpretation of this oscillationcan suggest an error larger than is true. For instance, if an optima were to exist approximately47mid-way between the two results, themargin of error between the upper and lower approximationwould be about double the actual deviation from the optima.3.2.3 Price and Load BlocksPrice and Load Blocks are the alternative approach to incorporating inflow-dependent marketprices and load projections without including an additional random state variable in the SDDPproblem formulation. Price and loadblocksprovidesmore accurate representation andoperationalpolicies regarding reservoir storage accounts. The creation of these price and loads blocks is a bitarduous, however they only need to be calculated once for any of the possible inflow years thatthe inflow model provides. For instance, if the possible seasonal volume forecasts are predictedto be similar to those of 1972, then the price and load blocks need only be generated from 1972 to1977, as the model is currently only being used to optimize the next five years of operation (or 60months). A flowchart demonstrating the procedure for generating price load blocks is presentedin Figure 3.6.First, a seasonal forecast is generated by the AR(1) inflow model and a matching water year isselected using the forecast seasonal volume river flows. For example, consider the year 1972. Next,General Operations Model (GOM) data provided by BC Hydro are pulled for the correspondingwater year. Since the SDDP model considers 60 stages, or five years, six years of total data arepulled. It is also important to remember that water years begin in October, not January. GOMdata provide optimal plant generation at hourly intervals for the optimized BC Hydro facilities.These GOM optimizations are again based off historic inflow sequences. An example of the GOMdata for 1972 is provided in Table 3.1. Note that some of the columns and rows were removedfor the sake of space. Once all six years of GOM load data have been pulled, the residual load iscalculated. The residual load is defined as the system load net of the optimized plant generation:R = BC_LOAD− (SEV_GEN+ALN_GEN+ARD_NT_GEN+NTS_GEN+ BC_WIND_GEN+ BC_IPP_GEN+BC_SHYDROGEN+ BC_THERMAL_GEN)(3.10)48Table 3.1: 1972 GOM Load Data (MW)t hours ARD_GEN GMS_GEN MCA_GEN PCN_GEN REV_GEN . . . BC_LOAD1 1 174.913 245.553 0 106.796 386.261 5582.912 1 174.89 245.554 0 106.796 163.759 5368.813 1 174.866 144.891 0 106.796 138.718 5246.134 1 174.842 95.97 0 106.796 138.719 5208.845 1 174.818 103.482 0 106.796 149.019 5218.466 1 174.799 245.557 0 106.796 411.792 5337.54...8760 1 176.51 476.132 1303.072 123.492 1149.634 . . . 6105.34The hourly residual loads are then separated into separate data structures, one for each of HLH(heavy load hours) and LLH (light load hours). HLH pertains to the residual loads between thehours of 06:00 to 21:00, inclusive, and LLH to 22:00 to 05:00. Next, these separated HLH andLLH data structures are used to sum the total HLH and LLH residual loads for each day (24 hourperiods). Note that t = 1 corresponds to midnight of the first day of the water year from the GOMdata. These summations are then used to generate residual load ratios by dividing the hourlyresidual loads by the corresponding daily total. For example, the first hourly residual load inOctober would be divided by the sum of the LLH (since midnight falls within the LLH coverage)for that day. These residual load ratios effectively serve as a means to see approximately whatproportion of the daily LLH or HLH are utilized at each hour.Next, daily load andprice predictions for 2021-2026 are retrieved fromavailable BCHydro data.BC Hydro has considered water years from 1972 to 2020 when developing predictions for a five-year window from 2021-2026. This consideration means that a unique 2021-2026 projection existsfor each inflow from all water years from 1972 to 2020. BC Hydro’s projection data is providedat a daily resolution and must be expanded to hourly if it is to be used in block calculations. Thecalculated price and load ratios (discussed below) are used to do this by multiplying each dailyprojection by a corresponding ratio for each given hour of the day. BC Hydro’s data is providedfor both LLH and HLH as separate datasets. As such, respective ratios are used when mappingdaily LLH and HLH predictions to hourly. This is discussed further in the following sections.49Figure 3.6: Price and Load Blocks Flowchart50SDDP Price BlocksPrice block calculations begin by first accessing market prices provided by the GOM data. Table3.2 shows a typical price file from the GOM data.Table 3.2: 1972 GOM Price Data (US$/MWh)t . . . US_IMP_PRICE_AVG US_EXP_PRICE_AVG . . . US_TIE_MIN_AVG1 83.95 75.504 -1880.92 83.95 75.504 -1880.93 80.042 71.989 -1880.94 77.755 69.932 -1880.95 83.95 75.504 -1880.96 90.897 81.753 -1880.9...8760 99.323 89.719 -1895.1Several columns have been removed for the sake of space. The removed columns include:US_Imp_SPOT_AVG, US_Exp_SPOT_AVG, and US_TIE_MAX_AVG. Prices used in block calcula-tions are determined by a wheeling calculation. A wheeling fee or price is an amount in $/MWhwhich a transmission owner recovers for the use of its system. This price is calculated simply by:W = US_IMP_PRICE_AVG−US_EXP_PRICE_AVG (3.11)These hourly GOM prices are sorted into HLH and LLH’s before their subsequent price ratiosare calculated. These ratios are calculated by dividing each hourly price by the average price forthat 24 hour period. These ratios are then multiplied by the HLH and LLH market predictionsprovided by BC Hydro to map the daily projections onto hourly. As mentioned previously, pricepredictions are based on water years from 1972-2020. Price blocks are broken down into fivecomponents, namely the peak price hours (PPH), highest HLH hours (HHLH), HLH hours, LLHhours, and lowest LLH hours (LLLH). These each correspond to the hours of greatest demand,second highest, third, second lowest, and lowest, respectively. The allocation of hours to eachdivision is typically 2-6-8-6-2, meaning PPH’s correspond to the highest two hours of demand,HHLH the next highest six, and so on. Each day of each month is sorted according to these51divisions, meaning each month (one group for each year, e.g. Dec 2021, Dec 2022, etc.) will havefive groups, one for each division. So, for each day in each month, the PPH, HHLH, HLH, LLH,and LLLH are computed and an average is taken. This average represents the price ’block’ for thatmonth. The 1972 price blocks are provided in Table 3.3 for reference.Table 3.3: 1972 Price Blocks (US$/MWh)BLOCK 1972 PricePeak HHLH HLH LLH LLLHDec-21 82.55 64.70 45.59 37.68 25.50Dec-22 89.47 70.50 56.88 50.47 29.22Dec-23 77.98 67.11 56.12 50.82 37.31Dec-24 75.87 64.95 53.77 47.22 27.65Dec-25 73.66 62.77 54.17 49.63 39.15Dec-26 102.07 79.71 66.56 60.63 46.50SDDP Load BlocksOnce hourly projections for all days from December 2021 to November of 2026 have been mappedfrom daily data, it is possible to calculate the load blocks. Load blocks are calculated identicallyto price blocks. An example of December load blocks for 2021 to 2026 from a 1972 inflow year isprovided in Table 3.4.Table 3.4: 1972 Load Blocks (MWh)BLOCK 1972 LoadPeak HHLH HLH LLH LLLHDec-21 5725.246 5444.617 4562.954 4336.313 3046.661Dec-22 4410.023 4171.562 3810.87 3529.93 3550.296Dec-23 4734.713 4759.288 4659.49 4437.76 4423.248Dec-24 4380.195 4266.549 4355.938 3564.543 3193.071Dec-25 5378.777 5210.954 5007.629 4317.984 5041.796Dec-26 5296.065 5201.473 4718.375 3630.455 4102.078It is interesting to note that lower demand hours do not always average out to be less than theirhigher counterpart. For example, in some instances in Table 3.4 HHLH is actually greater thanthe peak. The same can be seen for some of the LLLH blocks. The SDDP model uses blocks thatare sorted by price prior to determining the corresponding load, so this can explain the observed52behaviour in Table 3.4. It should be noted that henceforth, any use of the term ’load’ is referring toresidual load as it is used in the load blocks.3.2.4 Price and Load Block Implementation in the Julia LanguageExisting SDDP Julia ModelThe existing SDDP model considers only a single price and load block. This block corresponds tothe 1972 water year and was calculated in the same way discussed in §3.2.3. This means that forevery selected water year corresponding to the AR(1) generated inflow, the 1972 price and loadblocks were used regardless. This provides a very small parameter space for both price and loadwithin the model. The existing SDDP model implements this single year of blocks by referencingan external (to the primary model file) datafile. This datafile is referenced by the primary modelfile as needed at each iteration. Prices are included in the objective formulation as a component ofthe second term in (3.1). Load blocks serve to provide right-hand side limits on loading constraints.Forward, backward, and simulation passes always reference the 1972 price and load blocks, anddo not at any time update this datafile, regardless of the results of the inflowmodel. It is apparentprior to assessing results that the shape of the objective function as well as reservoir releases willbe tightly governed by the relatively limited range of prices and loads provided by one year ofblocks. Existing price and load block implementation is presented in Figure 3.7.Computational Tractability and the Master SheetFor any given year, using a personal computer to calculate the corresponding price and load blocksalone takes approximately 45 seconds. Of course, this will vary depending on the power of themachine used. A separate program from the SDDP model is used when calculating these blocks,and was originally going to be incorporated into the existing model to calculate blocks at everyiteration, dependent on the selected water year.53Figure 3.7: Current Price and Load Block Implementation with Model FlowHowever, consider an SDDP model run with 4000 iterations. If price and load blocks wererecalculated for every iteration, and if every calculation took 45 seconds, the model would take anadditional 125 days to run. This is not computationally tractable. For this reason, implementationinto the existing model was done using what can be referred to as a ‘master document’ or ‘mastersheet’. This master sheet was created by merging all price and load blocks calculated prior toimplementation. The master sheet contains price and load blocks from 1972 to 1999 as availableGOM data does not exist beyond 1999. The existing SDDP model algorithm requires 49 historic54sequence simulations (1972-2021) to check for convergence. This means that convergence checksfollowing the implementation of themaster sheetwill only utilize 27 (1972-1999). Thiswill certainlyhave effects on convergence trends, however it is not clear to what degree this reduction in historicsequence simulations would alter model results. A basic comparison was made by running theexisting model with both 27 and 49 year historic sequence simulations and comparing results.However, in an attempt facilitate more directly comparable results between the existing and newmodel implementation, the existing model convergence check was changed to also utilize only 27historic sequence simulations when assessing against the implemented changes.3.3 Existing Model Results3.3.1 ConvergenceThe results of the unchangedmodel are presented and assessedhere for both 49 years of simulation.The same model results for 27 years are presented in Chapter 4. This will provide insight intothe implications of reducing the simulation coverage, as well as existing model behaviour. Forboth renditions, 500, 1,000, 1,500, 2,000, 2,500, 3,000, and 4,000 SDDP iterations were conducted.Iterations beyond 4,000 caused substantial iteration slowdown with the implemented JuMP solver’Clp’. These issueswere related tomemory usagewithin the program. Rudimentary investigationssuggest this iteration slowdown to be related to the storage of model results in Julia memory untilall iterations were complete. Future investigations may wish to implement new code to ’dump’memory every 1500 iterations by writing results to a datasheet prior to full run completion. Figure3.8 shows the trend of the final simulation result for all considered iterations, i.e. the final objectivefunction value for each month for the 49th simulation.Objective value trends remained very close to each other for the first three years, diverging sig-nificantly in year four. Convergence was defined to be the difference between two final simulationresults of less than ε = 0.05 or 0.5%. This selection of criteria was made with knowledge of theproof of convergence from Linowsky and Philpott, 2005, wherein it is stated that after a sufficientlylarge number of cuts, an SDDP model which satisfies four criteria always has an optimal solution.This SDDP model satisfies all four criteria. These criteria are:55Figure 3.8: Final Simulation Objective Value Trends1. Random quantities appear only on the right-hand side of the linear constraints in each stage,2. The set of Ωt of random outcomes in each stage t is discrete and finite,Ωt = {wti|i = 1, . . . , qt < ∞} with probabilities pti > 0, ∀i,3. Random quantities in different stages are independent,4. The feasible region of the linear program in each stage is non-empty and bounded.Thus it follows that if all cuts are sufficiently close, they have approached the optimal solution. Fig-ure 3.8 demonstrates the previously mentioned oscillation of simulation results around an optima,specifically in year four when observing the 500, 3,000, and 4,000 iteration results. Approachingthis optima from two directions by producing Benders cuts which either overestimate or under-estimate the actual future value function. Provided enough iterations, results will still approachthe same optima assuming a feasible solution exists, which is the case here as the SDDP.jl packagerequires complete recourse and cannot continue when presented with an infeasible solution. It56should be noted that following analysis of the results presented in Chapter 4, convergence wasdetermined to not be a relevant metric when comparing the existing model. The overall modelresults were substantially different and thus it was more pertinent to compare storage, marginalwater value, and objective function results against each model only and not consider convergencetime against themselves. Figure 3.8 clearly demonstrates a poor implementation of the TerminalValue Function (TVF) in the existing model. The final three months result in the objective functionterminating at zero. The existing model TVF essentially values future generation outside the fiveyear planning horizon as zero. This results in a model which is attempting to ascertain massiveobjective function results without regard to the horizon following the five years. This is alsoapparent in the marginal water value and storage results presented in §3.3.2 and §3.3.3.The model came exceptionally close to convergence after 2,000 iterations, where the averageε was 0.23%. However, month 21 and 55 had an error of 0.83% and 0.7% respectively, renderingconvergence incomplete. Interestingly, subsequent iterations of 2,500, 3,000, and 4,000 provided anotable increase in error. Table 3.5 provides the average, maximum, and minimum errors betweensimulation results. Although all comparisons provided a minimum error of 0%, the monthswhich provided these did not necessarily always provide 0%. For example, An ε value of 0% wascalculated for the 36th month between the results of 2,000 and 1,500 iterations. This value thenjumped to 4.4% between 2,500 and 2,000 iterations. Note that the calculation of ε¯ did not takeinto account the final three months of simulation results, as all water values were zero, alwaysproviding an error of 0% which would provide a unrepresentative average if included. Table 3.6shows the error change for every month. The first row corresponds to the difference between 1,000and 500, second row for 1,500 and 1,000, etc. It is very clear that the convergence trends within thefourth year struggles.Table 3.5: ε Variations for 49 SimulationsIterations500 1,000 1,500 2,000 2,500 3,000 4,000ε¯ - 2.16% 0.84% 0.23% 1.20% 2.11% 2.37%Max ε - 10.23% 2.63% 0.87% 4.87% 13.74% 11.50%Min ε - 0% 0% 0% 0% 0% 0%57Table 3.6: ε by Month for 49 SimulationsDecember January February March April May June July August September October NovemberYear 10.20% 0.16% 0.15% 0.22% 0.60% 1.16% 0.92% 0.14% 1.13% 1.52% 1.53% 1.26%1.24% 0.94% 0.89% 0.48% 0.02% 0.62% 0.59% 0.52% 0.07% 0.34% 0.37% 0.56%0.29% 0.15% 0.14% 0.16% 0.17% 0.27% 0.87% 0.45% 0.22% 0.33% 0.26% 0.16%0.50% 1.55% 1.47% 1.37% 1.32% 1.14% 0.54% 0.25% 0.46% 0.36% 0.38% 0.26%0.16% 0.19% 0.19% 0.18% 0.17% 0.47% 0.61% 0.32% 0.13% 0.01% 0.00% 0.07%0.06% 0.67% 0.63% 0.60% 0.74% 1.29% 1.83% 1.06% 0.34% 0.24% 0.21% 0.29%Year 21.16% 1.08% 1.01% 0.97% 0.97% 1.70% 1.04% 1.11% 2.81% 2.04% 1.97% 1.76%0.58% 0.59% 0.60% 0.58% 0.47% 0.25% 0.31% 0.05% 0.79% 0.50% 0.45% 0.48%0.08% 0.04% 0.00% 0.01% 0.03% 0.15% 0.30% 0.36% 0.79% 0.40% 0.37% 0.36%0.19% 0.20% 0.20% 0.19% 0.25% 0.19% 0.48% 0.40% 1.08% 1.06% 1.03% 0.95%0.09% 0.08% 0.07% 0.07% 0.06% 0.07% 0.07% 0.03% 0.07% 0.46% 0.47% 0.46%0.27% 0.23% 0.22% 0.22% 0.15% 0.26% 0.14% 0.21% 0.20% 0.24% 0.25% 0.30%Year 31.59% 1.43% 1.35% 1.25% 1.30% 1.24% 1.17% 1.31% 1.53% 3.12% 10.23% 9.62%0.43% 0.37% 0.35% 0.33% 0.38% 0.04% 0.29% 0.45% 0.36% 0.29% 2.62% 2.63%0.34% 0.28% 0.26% 0.25% 0.08% 0.10% 0.17% 0.14% 0.17% 0.17% 0.17% 0.05%0.84% 0.78% 0.74% 0.70% 0.69% 0.87% 0.75% 1.06% 0.68% 2.08% 4.87% 4.40%0.42% 0.38% 0.36% 0.34% 0.18% 0.43% 0.03% 0.27% 0.22% 0.71% 0.57% 6.33%0.27% 0.25% 0.24% 0.23% 0.22% 0.52% 0.12% 0.05% 0.00% 0.02% 2.49% 8.64%Year 48.88% 8.33% 7.95% 7.57% 7.47% 7.99% 7.04% 0.64% 0.65% 0.69% 0.29% 0.29%2.44% 2.28% 2.16% 2.04% 1.98% 2.02% 2.14% 0.15% 0.78% 1.16% 1.29% 1.28%0.01% 0.04% 0.04% 0.04% 0.11% 0.05% 0.06% 0.18% 0.24% 0.35% 0.34% 0.33%4.05% 3.75% 3.56% 3.39% 3.43% 3.52% 3.75% 2.18% 0.29% 0.39% 0.41% 0.41%5.88% 5.54% 5.38% 5.20% 5.45% 6.61% 8.21% 11.13% 13.39% 13.74% 12.97% 8.44%7.97% 7.55% 7.28% 7.01% 7.16% 8.44% 9.93% 10.55% 11.42% 11.50% 11.27% 7.59%Year 50.20% 0.14% 0.03% 0.04% 0.03% 0.91% 0.38% 0.00% 0.00% 0.00% 0.00% 0.00%1.24% 1.10% 1.01% 0.92% 0.89% 0.59% 0.96% 0.00% 0.00% 0.00% 0.00% 0.00%0.29% 0.26% 0.27% 0.28% 0.25% 0.13% 0.67% 0.47% 0.00% 0.00% 0.00% 0.00%0.50% 0.51% 0.52% 0.54% 0.60% 0.68% 0.04% 0.47% 0.00% 0.00% 0.00% 0.00%0.16% 0.14% 0.05% 0.05% 0.09% 0.26% 0.09% 0.68% 1.32% 0.00% 0.00% 0.00%0.06% 0.13% 0.10% 0.13% 0.13% 0.01% 0.01% 0.76% 0.99% 0.00% 0.00% 0.00%583.3.2 Storage AccountsIt is also pertinent to consider the storage account operations of the reservoirs, specifically ARD,GMS, MCA, ARD treaty, and MCA treaty storage accounts. These plots were generated using aJulia script based on the output data of the model. These five plots are presented below. Each plotcontains the final simulation result same as previously mentioned.Figure 3.9: ARD Storage Trends (49 Sim Years)ARD storage trends follow similar behaviour to what was seen in the water value error trends.Storage accounts for each iteration tend to be extremely similar, if not identical. That is, until the 35thmonth, where the results of each iteration begin to change drastically. Until this point, all storageaccounts are exactly identical. During the 35th to 55th month, the differences in results range from11.15% to 276.77%. GMS appears to be immune to the factors influencing these deviations. GMSstorage accounts for every iteration are identical, with 0%deviations at every stage. Comparatively,MCA is behaves similar to AR in terms of deviations. As with ARD, until approximately the 35thmonth, results are identical. The same is seen for MCA treaty storage trends.59Figure 3.10: GMS Storage Trends (49 Sim Years)Figure 3.11: MCA Storage Trends (49 Sim Years)60Figure 3.12: ARD Treaty Storage Trends (49 Sim Years)Figure 3.13: MCA Treaty Storage Trends (49 Sim Years)61Figure 3.14: 1972 Price BlockFigure 3.15: 1972 Load BlockThese deviations do not appear to correlate to either the 1972 price or load block. Instead,it appears that these changes in water value are governed by storage policies surrounding theeffected months. This is very apparent when observing Figures 3.10 through 3.13. Figure 3.14 and62Figure 3.15 show the 1972 trends in the price and load blocks respectively. Figure 3.10 shows thesubstantial influence the TVF has on operating storage accounts at GMS. Asmentioned previously,since the TVF does not place any value on generation outside the five year planning horizon, theexisting model is attempting to obtain as much objective function value as possible prior to the60th month.3.3.3 Marginal Water ValueConvergence criterion was defined in terms of overall objective function trends rather than themarginal value obtained from Benders cuts (the slope of the hyperplane, βt,n) as no discernibletrend towards a converged future value functionwaspresent for eachmonth at each facility. Figures3.16 through 3.18 demonstrate the lack of a converging trend. It may seem at first glance that GMSand MCA have indeed converged, however, their differences between successive iterations remainalmost constant and did not improve noticeably from 500 to 4,000 iterations as shown in Table 3.7.Figure 3.16: ARDMarginal Water ValueWhen considering the subsequent Benders cuts resulting from the slopes provided in Figure633.16, the discrepancies between the 3,000 and 4,000 iterations far exceeded any reasonable conver-gence criterion for error. When assessing the cuts determined for GMS and MCA, they do indeedbehave much better. However, to consider this as a means for convergence would require that allfacilities were well behaved, which is simply not the case for ARD. Even then, regardless of thebehaviour or similarities between marginal water values at any given iteration, Figures 3.17 and3.18 once again provided trends which are essentially meaningless due to the TVF. The marginalvalue of water quickly reaches zero in the final three months as once again, the existing model isattempting to extract as much value as it can prior to the end of the five year planning horizon.This lack of consideration past the five year horizon is assessed in more detail in Chapter4. Again,the apparent convergence in marginal water values is essentially meaningless given the apparentover-utilization of storage accounts.Figure 3.17: GMS Marginal Water Value64Figure 3.18: MCAMarginal Water ValueTable 3.7: GMS and MCAMarginal Value TrendsIterations500 1000 1500 2000 2500 3000 4000GMSε¯ - 0.06% 0.08% 0.07% 0.45% 0.50% 0.22%Max ε - 3.40% 3.29% 3.40% 3.29% 25.76% 0%Min ε - 0% 0% 0% 0% 0% 0%MCAε¯ - 3.77% 2.64% 2.07% 2.87% 4.13% 5.30%Max ε - 36.25% 18.30% 12.72% 13.36% 15.28% 33.61%Min ε - 0% 0% 0% 0% 0% 0.02%653.3.4 Implications of 27 vs 49 Simulation YearsThe behaviour of the existing model with only 27 years of simulation was surprisingly differentthan that of 49. Storage operations had not yet approached similar results after only 27 years.Figures 3.20 through 3.21 demonstrate that for all storage accounts other than GMS, the simulationresults had not yet reached those presented in §3.3.2.Figure 3.19: GMS Storage Trends (27 Sim Years)66Figure 3.20: ARD Storage Trends (27 Sim Years)Figure 3.21: MCA Storage Trends (27 Sim Years)67Figure 3.22: ARD Treaty Storage Trends (27 Sim Years)Figure 3.23: MCA Treaty Storage Trends (27 Sim Years)68This leaves room for the question of how to improve the quantity of simulation runs while stillmaintaining the inflowdependencyof eachprice and loadblock. It is possible to extend simulationsindefinitely, however the simulation pass would need to become disconnected from the price loadblocks to do so. This serves to show that the variance in the new model results presented in§4.3 with respect to storage accounts is likely not only due to the increased parameter-spaces, butpossibly in part to the reduction in simulation passes.3.3.5 Computational EffortThe required computation time to complete model runs with linearly increasing iteration require-ments does not scale 1:1, but instead scales proportionately to the power function y = axb. Figure3.24 shows this relationship between computational time and iterations. A secondary plot ’rela-tive’ is provided within the figure which shows the time requirements relative to the time for 500iterations, as the precise time requirement in seconds will not be the same across all machines. Forexample, from Figure 3.24 1000 SDDP iterations takes approximately 250% longer than 500 SDDPiterations. The line of best fit which projects the expected computational effort per iteration is thefunction:Figure 3.24: Computation Time for Iteration Requirements69y = 0.0052x1.8961x (3.12)where y represents the time in seconds and x the number of iterations. This line of best fitprovides an R2 value of approximately 0.97 However, as additional data points are added, theaccuracy of (3.12) is reduced. Based on the available data points, previously representative powerfunctions when considering one less data point proved to always be an underestimate of the actualtime requirement. This projection is provided in Figure 3.25. According to (3.12), it would takeapproximately 55 hours to achieve 10,000 iterations. However, personal experience at attemptsof 10,000 iterations suggest a time considerably longer, as crashes would occur at approximately5,500 iterations and an elapsed time of 20 hours. It is not believed that these crashes were relatedto model infeasibilities, but rather memory limitations within Julia. Previously, infeasibilities hadbeen encountered and were handled by JuMP’s built-in error handling rather than crashing.Figure 3.25: Computational Time Projections70Chapter 4Results and Discussion4.1 Implemented ChangesWater year (or inflow) dependent price and load blocks were implemented through the use of a’Master Sheet’. This document contained all data regarding price and load projections with respectto the water years 1972-1999. The use of a master document was decided due to the computationalrequirements of generating new price and load blocks dynamically each time they are required.The time to generate a single set of blocks ranges anywhere from 45 to 60 seconds, dependent onthe computer used.Additional functions were added to the forward, backward, and simulation passes withinthe model. These functions performed only two tasks: limiting the water year selection to bebetween 1972 and 1999, and creating a new temporary data file for the corresponding price andload blocks from the master document. The new model flow is shown in Figure 4.1. Whileeffective, this implementation is not entirely efficient. The current structure of the Julia files limitsthe scope of both the constraint and objective function models to their respective files. This meansthe the master document cannot be directly referenced with respect to the generated water yearfrom the constraints or objective function files as they do not have scope access to the forward,backward, or simulation passes. Scope in computer science refers to an area of a program wherean item that has an identifier name is recognized. Unlike Python, Julia cannot access specificfunctions from other files by directly importing them unless the function is part of a module. Theversion of the SDDP model used here does not yet incorporate modules. Future improvementscan be made to this implementation should refactoring occur. Refactoring refers to the process71of restructuring existing computer code without changing its external behavior. As this modelis undergoing constant development internally at BC Hydro, it is very likely that by the time theresults presented here are available, they are out of date.Figure 4.1: Model After Implemented ChangesInitial implementations of this new model flow required over an 800% increase in computationaleffort per iteration, placing 4,000 iterations at over 80 hours to achieve, however the final imple-mentation required only an initial 76% increase, decreasing to less than 10% as the number ofiterations increased. This is discussed further in §5.1.72Figure 4.2: New Price Parameter-Space vs. 1972Figure 4.3: New Load Parameter-Space vs. 1972The parameter-space for market price and load demandswas substantially increased as a resultof the master document use. The ranges of prices and loads expanded by up to 123.8% and 31.3%respectively. The maximum and minimum for each block component was determined for each73month across all price and load blocks. These points were used to compare the new maximumand minimum parameter-spaces. The 1972 price and load blocks are shown alongside as a pointof reference. Figures 4.2 and 4.3 show the spread of prices with respect to peak hours for price andload blocks. The unshaped GOM data associated with each block is also shown. Expansions forHHLH, HLH, LLH, LLLH were omitted for the sake of space.Additionally, the terminal value function was updated such that the end of planning horizonwater value does not terminate to zero. Instead, the results of several years of simulationwere usedto determine an appropriate TVF. These simulations and the final selection of TVFwere completedby Dr. Moosavian of BC Hydro. As such, the final results presented in this Chapter will not onlybe influenced by the implementation of price and load blocks, but also the addition of a moreappropriate TVF.4.1.1 Challenges EncounteredFigure 4.4: Negative Load in LLH BlocksOccasionally negative values were calculated during the creation of a price or load block compo-nent. This was much more common with price blocks than load. To rectify this, the negative loadand price requirements were instead defaulted to zero. This solved a bounding issue experienced74by the Clp solver as it would misinterpret negative loading as a non-existent lower bound. Similarbehaviour was observed at blocks with negative price. The load space which provided the mostsignificant negative results was observed to be the LLH blocks as seen in Figure 4.4.Reiterating what was mentioned in §4.1 with respect to code factorization, the implementationof thenewmodel can certainly be improved. However, the additional computational effort requiredper iterationwas not increased in such away thatwould render themodel impractical to run. Thus,therewas no attempt to improve the efficiency of the entiremodel. Efficiencieswere instead focusedon the additional three functions which were implemented.4.2 Objective Function TrendsThe results of the new model when compared to the existing model with new TVF followed asubstantially different trend than those presented in the from Figure 3.8. The most apparentchange in the objective function results is the seasonally-consistent final water values. Figure 4.5plots both the results from the existing model at 4,000 SDDP iterations and the new results at 5,000SDDP iterations. Note that all results in this section which contain dollar values in their units arenormalized for confidentiality. Comparisons are thus made only based on percentages.As expected, objective function results from both models provide seasonal variations in watervalue. Existing model results provided surprisingly similar results for all months except thosebetween month 41 and 51. Within these ten months, the exiting model results present andapproximately 65.5% decrease in water value. The overall average difference between results ateach stage is 69.6%. Lowest values for bothmodels occur aroundmonth 51, before quickly peakingprior to month 60. This change from the existing model without an appropriate TVF shows itsimpact on final planning horizon results. Overall, the total water value extracted by the existingmodel with the new TVF is greater than the newmodel by approximately 7%. This is likely due tothe stricter market data used to generate the Bender’s cuts at each stage within the existing model,resulting in an over-estimation. A 7% differencemay seem small, however due to themagnitude offinances involved, this difference has significant potential to impact operational decisions. Whilethis comparison between final objective function results can be informative, it must be rememberedthat drawing a conclusion from only the comparison of 1972’s data will not fully encompass the75Figure 4.5: NewModel Final Simulation Objective Value Trendsexisting models behaviour. It is entirely possible that if data for any other year was used, theexisting model may have instead under-estimated water value when compared to the new model.4.2.1 ConvergenceWhile consideration of convergence is important, a comparison between the existing and newmodels was not conducted. The implementation of the TVF brought with it significant changesthroughout all results. As such, a comparison was considered more or less irrelevant at this stage.It is likely that model parameters (such as the TVF and price and load blocks in this case) shouldbe assessed for overall impact prior to assigning much merit to convergence. In practice, SDDPmodels are often run for a smaller number of simulations such that only the first year or two ofresults are sufficiently converged. This is of course context dependent, however it is important tonote that convergence at each stage is not necessarily required for aspects of themodel to be usable.764.3 Storage AccountsFigures 4.6 through 4.10 present the storage account results from the new model compared toexisting. GMS appears to follow very similar operations during the five year planning horizon.Noticeable differences are present between months 18 to 50, where the existing model utilizesmuch less of the available storage account. The largest deviation occurs at month 22, where theexisting model presents a 26% decreased GMS storage utilization. A significant drop in storageis observed between months 40 and 50, much in line with the reduction in objective functiontrends around the same time. It is unlikely that this reduction in storage utilization is causing thereduction in objective function results, as even though this trend is present in both the existingand newmodels, it is only the existing model that sees a significant objective function drop withinthese months. A sensitivity analysis regarding specifically GMS utilization and its effects on watervalue could be conducted to better understand this correlation.Figure 4.6: NewModel GMS Storage TrendsOpposite to the very seasonal variations observed for GMS, ARD storage trends aremuchmoreconservative, spending many months at maximum capacity. This is true for both the existing and77new model. Overall, no significant or noteworthy difference are observed.Figure 4.7: NewModel ARD Storage TrendsNew MCA storage accounts are very in seasonality to the existing model but differ greatly inmagnitude, especially during the final 30 months. Additionally, peak values appear to be offset byapproximately three months, with the existing model providing peaks following the new model.This gap appears to also grow larger at the later months. Between months 30 and 50, existingstorage accounts are up to 48.5% larger. These differences are possibly due to the variances instorage limitations imposed by the load blocks.TY storage accounts for bothARDandMCAshowvery similar results frombothmodels. Recallthat negative storage accounts observed in Figure 4.10 are possible due to the implementation ofpenalty functions. A negative value here does not reflect a physical storage volume, but ratherrepresents values below a minimum threshold.78Figure 4.8: NewModel MCA Storage TrendsFigure 4.9: NewModel ARD Treaty Storage Trends79Figure 4.10: NewModel MCA Treaty Storage Trends4.4 Marginal Water ValueMarginal water values are compared to assess the overall trend of the Benders cuts at each stage.Recall themarginalwater values in thiswork are calculated as the slopes of thewater value function,and are used to approximate the future value functions at each stage. Results of this section arevalidated using a Successive Linear Programming (SLP)model of the same system. This validationis presented in §4.5.Figures 4.11 through 4.15 present the trend of marginal water values for each facility. GMSmarginal value trends demonstrate a drastically more consistent trend from the new model. Thejagged peak pattern presented by the existing model without a TVF (as seen in Figure 3.17) areno longer seen within the new model. Additionally, the concerning ’termination to zero’ trendobserved in the existingmodel beginning atmonth 57 is no longer present. The largest discrepancyfor GMS occurs at month 41, where the existing model results present a marginal water value ofzero. Again, this could be a component in the cause of the low water value presented duringsimilar months in Figure 4.5.80Figure 4.11: NewModel GMS Marginal Water ValueFigure 4.12: NewModel ARD Marginal Water Value81Figure 4.13: NewModel MCAMarginal Water ValueFigure 4.14: NewModel TY ARDMarginal Water Value82Figure 4.15: NewModel TY MCAMarginal Water ValueARDmarginalwater is generally consistent betweenbothmodels at all stages other thanmonths44 to 50. Within these months, the existing model presents a marginal water value which is a 70%decrease when compared to the new model. As well, the final three months of the existing modelagain do not terminate to zero, and continue to behave as would be expected seasonally.MCA trends stay consistent which what has been observed thus far. Both the existing and newmodel provide relatively similar results, with the largest discrepancy occurring at month 47.Marginal water values for both MCA and ARD do not provide much in terms of any usefultrend. There does however seem to be a more seasonal trend in the spikes of the marginal watervalues. Both treaty accounts for both the existing and new model provide higher marginal watervalues from the new model. Treaty accounts marginal value peaks increased by approximately44% and 70% for ARD and MCA respectively.Overall, marginal value trends do not deviate by a significant amount. This could indicate thatbothmodels are developing similar future value functions shapes, albeitwithdifferentmagnitudes.This is an indication that the implementation of price and load blocks, while altering the coverageof the model, does not noticeably alter or hinder the development of future value functions.834.5 Validation of ResultsResults of the new model were checked alongside those provided by a successive linear programapplied to the same optimization problem. Successive linear programming is an optimizationtechnique that involves breaking a complex problem down into smaller, more manageable sub-problems, each of which can be solved using linear programming. This allows for an efficient anditerative approach to solvingdifficult optimization problemswhichmaybe nonlinear. In successivelinear programming, each sub-problem is solved using the solution from the previous sub-problemas the startingpoint, which canhelp to improve the overall accuracy and convergence of the solution(Nocedal & Wright, 2006). This technique is often used in a variety of fields, including operationsresearch, economics, and engineering, to solve complex optimization problems.Successive linear programming and stochastic dual dynamic programming are both optimiza-tion techniques that can be used to solve complex problems. However, there are some key differ-ences between the two.Successive linear programming involves breaking a complex problem down into smaller sub-problems, each of which is then linearized and solved using linear programming. This allows foran efficient and iterative approach to solving difficult optimization problems. In contrast, stochasticdual dynamic programming is a more general optimization technique that can be used to solve awide range of optimization problems, including both linear and nonlinear problems.Another key difference between the two techniques is the way in which they approach theproblem. In successive linear programming, the solution from each sub-problem is used as thestarting point for the next sub-problem, which can help to improve the overall accuracy andconvergence of the solution. In contrast, stochastic dual dynamic programming uses a dynamicprogramming approach, in which the solution is built up incrementally by making a series ofdecisions at each stage of the problem.Overall, while both techniques can be used to solve complex optimization problems, they differin the way they approach the problem and the types of problems they can handle. Successivelinear programming is specifically designed for solving problems using linear programming,while stochastic dual dynamic programming is a more general optimization technique. However,the results of both the SLP and SDDP model applied to the same system should provide relatively84similar results and can be used as a means to validate the development and progression of thenew SDDP model. The SLP results are presented from a 5,000 iteration model maintained by BCHydro.4.5.1 Storage AccountsFigures 4.16 to 4.20 compare the results of the new SDDP model to the established SLP model.The trends between both models for GMS storage remain similar. Peaks align within one to twomonths, with maximum and minimums deviating by no more than 24.4% and 23.6% respectively.ARD storage values also adhere to similar values from both models. The new SDDP modelresults appear to behave as an ’upper bound’ on those provided from SLP. SLP results are muchmore varied month by month. Recall however that the comparison between SLP and SDDP is notone-to-one, but rather is done only to provide insight into the progress of model development andvalidity of results. Figures 4.17 (ARD) and 4.19 (TY ARD) follow the same trends as do the SLPresults.Figure 4.16: SLP vs New SDDP GMS Storage85Figure 4.17: SLP vs New SDDP ARD StorageFigure 4.18: SLP vs New SDDP MCA Storage86Figure 4.19: SLP vs New SDDP TY ARD StorageFigure 4.20: SLP vs New SDDP TY MCA Storage87MCA storage accounts are larger from the SLP model than the new SDDP. Important trendsare similar between both results, namely the timing of peaks and toughs. However SLP does notpresent a monthly reduction in storage utilization similar to the new SDDP results. However,Figure 4.20 demonstrates a much closer relationship between SLP and SDDP results for MCAtreaty storage.Overall, the results provided by the new SDDP implementation behave well when comparedto SLP results. This improvement is even more noticeable if comparisons are made between theexisting model without the new TVF and SLP (not included here). It is important to reiterate that,while these two optimization models are being used to compare results, deviations are absolutelyexpected. Both models approach the system differently and make different assumptions, considerdifferent parameters, and implement different levels of randomness.4.5.2 Marginal Water ValueFigures 4.21 to 4.25 compare the results of the new SDDP model to the established SLP model.Figure 4.21: SLP vs New SDDP GMS Marginal Value88Figure 4.22: SLP vs New SDDP ARDMarginal ValueFigure 4.23: SLP vs New SDDP MCAMarginal Value89Figure 4.24: SLP vs New SDDP TY ARDMarginal ValueFigure 4.25: SLP vs New SDDP TY MCAMarginal Value90Discussion pertaining to each Figure is not required here as the only metric of interest is themagnitude of that marginal water value between both models. A direct comparison of trends isnot relevant as each model calculates marginal value much differently. Overall, for all facilitieswith the exception of MCA’s treaty account, the magnitude of the marginal water values betweenboth models remains adequately similar. Concerns had previously arisen regarding the marginalwater value pertaining to GMS and MCA, specifically that they seemed too large. However, thecomparisons presented in Figures 4.21 and 4.23 clearly show that the new SDDP implementationproduces results which are of an appropriate magnitude. MCA treaty marginal water values varysignificantly between SLP and SDDP models as seen in Figure 4.25.4.6 Computational EffortPrior to code optimization, newSDDP iterations took approximately 800% longer to completewhencompared with the existing model. The Julia code was subsequently optimized to allow quickeraccess to the price and load block master sheet at each iteration by ensuring all model componentscould access the data directly. However, existing model code should still be refactored to facilitatea simpler set of scopes. As previously mentioned, this could be achieved through implementationof Julia modules.The implemented code optimization reduced the additional computational requirement from800% to only 10%. For SDDP iterations below500, the new implementation requires 80%additionalprocessing time. However, for SDDP iterations at or above 4,000, the additional processing time isless than 10%. This behaviour is presented in Figure 4.26.91Figure 4.26: New and Old Computation Time by Iteration92Chapter 5Conclusions and Recommendations5.1 Discussion of ResultsThe current implementation of a single water year to determine price and load blocks does notprovide a sufficient parameter-space to generate water values which are reflective of reality. Ad-ditionally, the existing TVF does not account for value associated with generation for months oryears following the end of the planning horizon. This results in the existing model frequentlyreleasing all storage accounts and subsequent marginal water values terminating to zero duringthe final three months of operation. This is likely in part the cause of the consistent over-valuingof available water year over year. Thus, a TVF was added to the existing model alongside theimplementation of yearly price and load blocks. Recall that the ’existing model’ referenced here isindeed referring to the iteration which includes the TVF.Objective function trends provided by the new SDDPmodel show seasonal peaks whichmatchthose determined by the existing model. The existing model without a TVF provides peaks whichare offset by approximately three months. The existing model water value for the first peak atmonth five was 88% of the maximum observed water value (normalized). First peak water valuesfor the new model were 77%, a difference of 11%. Lowest values at troughs for the existing modeloccur at approximately month 51 with 19%. Lowest values at troughs for the new model occurat approximately month 51 as well with 5%, a difference of 14%. On average, the existing modelprovided water values 7% greater than the new model.Storage accounts were utilized much more actively as a result of the implemented changes.GMS responded more actively to seasonal changes. ARD utilization maximums were decreased93such that storage volume never reached zero. MCA reservoir utilization maximum decreased byup to 52.2%. Treaty storage for ARD remained closely linked to the operations of ARD. Treatystorage forMCAproduced less variance but did trend towards preferring to operate outside of CRTregulations by slightly increasing the duration and quantity of volume levels below the minimumamount.Marginal water value trends for all facilities remained relatively similar. GMS trends behavedwell between both models, with the existing model being slightly more volatile and reaching avalue of zero between months 41 and 44. ARD marginal water values peaks were increased byup to 100%. A more clear seasonal variance is observed within ARD. MCAmarginal water valuesdecreased for a majority of months. The magnitude of peaks decreased by approximately 26.2%.Treaty ARD and MCA marginal water value peaks increased by approximately 44% and 70%respectively. Treaty marginal water values appear to have developed a more seasonal trend.5.2 Recommendations for Future workFurther investigation into the development and implementation of a dynamic TVF could greatlybenefit future model results. The implications of assigning value to indefinite months of futuregeneration may provide results further reflective of reality. Additionally, continued improvementand refactorization of existing Julia code may result in faster model run times. Investigation intothe possible use of an alternative solver to JuMPmay be warranted. A detailed sensitivity analysismay provide useful information into the direct contributions each facility plays in model outputs.Based on the results presented here, it seems that GMS and MCA operations provide the largestimpactwith respect to objective function values. Overall, the consideration of price and load blocksby the new SDDP model produces results which seem more realistic, and future renditions of theSDDP model within BC Hydro should continue to consider their implementation.94BibliographyBenders, J. F. (1962). Partitioning procedures for solving mixed-variables programming problems.Numerische Mathematik, 4(1), 238–252. https://doi.org/10.1007/bf01386316Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge University Press, 216. https://doi.org/10.1017/cbo9780511804441Brewka, G. (1996). Artificial intelligence—a modern approach by stuart russell and peter norvig,prentice hall. series in artificial intelligence, englewood cliffs, nj. The Knowledge EngineeringReview, 11(1), 78–79. https://doi.org/10.1017/s0269888900007724British Columbia, G. o. (1964). The columbia river treaty and protocol.Department of External Affairs/ Northern Affairs and National Resources.C.E.R. (2018). Canada’s energy future 2018: Energy supply anddemandprojections to 2040.NationalEnergy Board. https://doi.org/https://doi.org/10.35002/wp48-nq03Chen, Z. L., & Powell, W. B. (1999). Convergent cutting-plane and partial-sampling algorithmfor multistage stochastic linear programs with recourse. Journal of Optimization Theory andApplications, 102(3), 497–524. https://doi.org/10.1023/a:1022641805263Conejo, A. J., Castillo, E., Minguez, R., & Garcia-Bertrand, R. (2011). Decomposition techniques inmathematical programming: Engineering and science applications. Springer.Dantzig, G. (1951). Application of the simplexmethod to a transportation problem.ActivityAnalysisof Production and Allocation, 359–373.Donohue, C. (1996). Sotchastic network programming and the dynamic vehicle allocation problem.Euclid. (1996). Euclid’s elements (Vol. 1). Department of Mathematics; Computer Science.Ferguson, T. S. (1989). Who solved the secretary problem? Statistical Science, 4(3). https://doi.org/10.1214/ss/117701249395Firsova, I. A., Vasbieva, D. G., Litvinov, A. V., Chernova, O. E., & Telezhko, I. V. (2019). Trends inthe development of the global energy market. International Journal of Energy Economics andPolicy, 9(3), 59–65. https://doi.org/10.32479/ĳeep.7728Fullner, C., & Rebennack, S. (2021). Stochastic dual dynamic programming and its variants. Karl-sruhe Institude of Technology.Guan, Z., Shawwash, Z., & Abdalla, A. (2018). Using sddp to develop water-value functions fora multireservoir system with international treaties. Journal of Water Resources Planning andManagement, 144(2). https://doi.org/10.1061/(asce)wr.1943-5452.0000858Iliadis, N., Tilmant, A., Chabar, R., Granville, S., Pereira, M., Consultoria, P., & Patria, R. (2008).Optimal operation of hydro-dominated multireservoir systems in deregulated electricitymarkets.Kazempour, J. (2022). Advanced optimization and game theory for energy systems.Kennington, J., & Helgason, R. (1984). Algorithms for network programming.Lillis, K. (2014). The columbia river basin provides more than 40% of total us hydroelectric gener-ation.Linowsky, K., & Philpott, A. B. (2005). On the convergence of sampling-based decomposition algo-rithms for multistage stochastic programs. Journal of Optimization Theory and Applications,125(2), 349–366. https://doi.org/10.1007/s10957-004-1842-zLucas, J., Sun, S., Zemel, R., & Grosse, R. (2019). Aggregated momentum: Stability through passivedamping. International Conference on Learning Representations. https ://openreview.net/forum?id=Syxt5oC5YQMurphy, J. (2013). Benders, nested benders and stochastic programming: An intuitive introduction.Cambridge University Engineering Department Technical Report. https://doi.org/10.48550/arXiv.1312.3158N.C.E.I. (2022). Billion-dollar weather and climate disasters. https://www.ncei.noaa.gov/access/billions/Nocedal, J., & Wright, S. J. (2006). Numerical optimization (2nd ed.). Springer-Verlag.Pereira,M.V., & Pinto, L.M. (1991).Multi-stage stochastic optimization applied to energy planning.Mathematical Programming, 52(1-3), 359–375. https://doi.org/10.1007/bf0158289596Philpott, A., & Guan, Z. (2008). On the convergence of stochastic dual dynamic programming andrelated methods. Operations Research Letters, 36(4), 450–455. https://doi.org/10.1016/j.orl.2008.01.013Rahmstorf, S., Foster, G., & Cahill, N. (2017). Global temperature evolution: Recent trends andsome pitfalls. Environmental Research Letters, 12(5), 054001. https://doi.org/10.1088/1748-9326/aa6825Rawas, F. (2020). Competition vs cooperation: Application of game theory in the multi-agentcoordination of a bc hydropower system. Faculty of Graduate and Postdoctoral Studies.Rougé, C., & Tilmant, A. (2016). Using stochastic dual dynamic programming in problems withmultiple near-optimal solutions. Water Resources Research, 52(5), 4151–4163. https://doi.org/10.1002/2016wr018608Slater, M. (2013). Lagrange multipliers revisited. Traces and Emergence of Nonlinear Programming,293–306. https://doi.org/10.1007/978-3-0348-0439-4_14Tret’yakov, A. A. (2009). Elementary proof of classical necessary optimality conditions for a math-ematical programming problem with inequality constraints. Doklady Mathematics, 79(2),243–245. https://doi.org/10.1134/s1064562409020252Vanderbei, R. J. (2004). Linear programming: Foundations and extensions. Kluwer.Wada, Y., van Beek, L. P., & Bierkens, M. F. (2011). Modelling global water stress of the recent past:On the relative importance of trends in water demand and climate variability. Hydrologyand Earth System Sciences, 15(12), 3785–3808. https://doi.org/10.5194/hess-15-3785-2011
Thesis/Dissertation
eng
Civil Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Shawwash, Ziad K. Elias
Parameter-space dependence in a stochastic dual dynamic programming model of a multireservoir hydropower system
Text