t l THE LONG-TERM OPERATION OF OVERYEAR STORAGE HYDROELECTRIC PROJECTS WITHIN A MIXED HYDROTHERMAL GENERATING SYSTEM by WILLIAM F. CASELTON B.Sc, University of Leeds, 1957 M.A.Sc, University of British Columbia, 1971 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of Civil Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH September, 1973 COLUMBIA In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the Head of my Department or by h i s representatives. It i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of The University of B r i t i s h Columbia Vancouver 8 , Canada Date ABSTRACT It is necessary to establish a new operating policy when a major hydro electric project is introduced into an existing generating system. When the project has a large reservoir, capable of storing substantial volumes of water over a number of years, then the operating policy cannot be confined to the current years operation - a long term operating policy is required. Analysis of the operation of the hydro project to determine the optimum long term policy becomes complex when the benefits of complementary interaction of the project with the rest of the generating system are considered. Further complication arises in the analysis from the uncertainty assoc-iated with reservoir inflows over this long time period. These complications prohibit the determination of an optimal long term operating policy by direct application of presently available optimisation techniques. In this thesis a system decomposition approach is presented for the determination of efficient long term operating policies for one or more overyear storage hydro projects in a mixed hydro and thermal electric generating system. The system is fir s t decomposed into a set of subsystems which are designed to preserve the complimentary operation of individual hydro and thermal plants. The long term operation i i i i i of a subsystem is optimised over a range of annual firm energy outputs using an analytical procedure which simultaneously applies optimisation techniques to both short term and long term operation. When only one overyear storage project exists in a generating system the results obtained from the analysis of its associated subsystem provide the required long term operating- policy. In the case where a number of overyear storage projects exist, or are proposed, i t is necessary to fi r s t analyse the associated subsystems individually. The results obtained form the basis of a procedure for determining an efficient operating policy which, over the long term, coordinates the operation of the overyear storage projects and at the same time integrates their operation with thermal and short term storage plants within the generating system. TABLE OF CONTENTS Page LIST OF TABLES LIST OF FIGURES . . . . . . . . . . . . CHAPTER 1. INTRODUCTION AND SUMMARY 1 1.1 GENERAL . 1 1.2 OPERATION OF HYDROELECTRIC SYSTEMS . . 2 2. DECOMPOSITION OF THE GENERATING SYSTEM. . . . . . . . . . . . . . 8 2.1 OPTIMIZATION OF COMPLEX SYSTEMS 8 2.2 DECOMPOSITION OF A MIXED HYDRO AND THERMAL GENERATING SYSTEM IN THE FORM OF A HIERARCHICAL SYSTEM. . . . . . . . . 13 2.3 SUBOPTIMIZATION OF OVERYEAR STORAGE SUBSYSTEMS . . 19 2.4 CO-ORDINATION OF OVERYEAR STORAGE SUBSYSTEMS 20 3. LONG-TERM OPERATION OF THE HYDROTHERMAL SUBSYSTEM . . . . . . . . 21 3.1 INTEGRATED OPERATION OF THERMAL AND HYDRO GENERATING PLANTS 21 3.2 INFLOW PREDICTION FOR THE OVERYEAR STORAGE HYDRO PROJECT . . 24 3.3 THE INDEFINITE TIME HORIZON 25 3.4 ASSESSMENT OF LOAD RELIEF DECISIONS UNDER UNCERTAINTY. . . . 27 3.5 SUBOPTIMIZATION SCHEME FOR THE HYDROTHERMAL SUBSYSTEM. . . . 29 4. ANALYSIS OF ONE YEAR OPERATION OF THE SUBSYSTEM. 30 4.1 DYNAMIC PROGRAMMING FORMULATION . 30 4.2 RECOVERY OF OPTIMAL STATE TRAJECTORY AND ASSOCIATED COST . . 34 5. FORMULATION OF THE ANNUAL LONG-TERM OPERATION OF THE HYDRO THERMAL SUBSYSTEM . . . . . . . . . . . 36 5.1 MARKOV PROCESS (OR CHAIN) IN DISCRETE TIME . . . . . . . . . 36 5.2 MARKOV PROCESS WITH RETURNS. . . . . . . . . . . . 38 5.3 MARKOV DECISION PROCESS 38 5.4 DECISION PROCESS WITH RANDOM INPUT 39 5.5 LONG-TERM OPERATION OF THE HYDROTHERMAL SUBSYSTEM AS A MARKOV DECISION PROCESS 43 5.6 DETERMINATION OF THE OPTIMAL LONG-TERM OPERATION POLICY FOR THE HYDROTHERMAL SUBSYSTEM . 44 iv V CHAPTER Page 6. COMPUTER PROGRAM AND COMPUTATIONAL EXPERIENCE. . 50 6.1 MAIN PROGRAM 50 6.2 DYNAMIC PROGRAMMING SUBROUTINE. . . . . . . . . . . . . . . . 53 6.3 SIMULATION SUBROUTINE . . . . . . . . 56 6.4 COMPUTATIONAL EXPERIENCE 58 7. STATISTICAL DATA FOR THE PORTAGE MOUNTAIN SUBSYSTEM 62 7.1 TEST FOR SERIAL DEPENDENCE OF ANNUAL INFLOWS. . . . . . . . . 62 7.2 PROBABILITY DISTRIBUTION OF THE ANNUAL INFLOW VOLUME. . . . . 63 7.3 DETERMINATION OF MONTHLY INFLOWS FROM A TOTAL ANNUAL INFLOW . 68 7.4 MONTHLY LOAD DISTRIBUTION 70 8. RESULTS OF THE ANALYSIS OF THE LONG-TERM OPERATION OF THE OVERYEAR STORAGE SUBSYSTEM. . . . . . . . . . 72 8.1 ANALYSIS OF THE HYDROTHERMAL SUBSYSTEM WITHOUT THERMAL CAPACITY CONSTRAINT . . . . . . . 72 8.2 ANALYSIS OF THE HYDROTHERMAL SUBSYSTEM WITH THERMAL GENERATING CONSTRAINT . . . . . . . . . . . . . . 82 9. CO-ORDINATION OF THE LONG-TERM OPERATION OF THE HYDROTHERMAL SUBSYSTEMS AND SUBGROUPS 89 9.1 CO-ORDINATION OF OVERYEAR STORAGE SUBSYSTEMS 89 9.2 CO-ORDINATION OF OVERYEAR STORAGE AND SEASONAL SUBGROUPS. . . 96 9.3 RECOVERY OF THE OVERYEAR STORAGE HYDRO PROJECT OPERATING POLICIES 99 9.4 CONCLUDING SUMMARY. . . . . . . . . . . . . . 1 0 0 REFERENCES • o * o o • « ' * o o e o o o o e o o o o o o • e o o o a o o o 102 APPENDICES o a » o « o o « o o o o o 9 o o o o o O o o o o o o s e o o « 1 0 4 LIST OF TABLES TABLE Page 7.2.1 RESULTS OF FITTING PEARSON TYPE I CURVE TO PORTAGE MOUNTAIN HISTORIC ANNUAL INFLOW VOLUMES 66 7.2.2 DISCRETE ANNUAL INFLOW VOLUMES AND PROBABILITIES 67 7.3.1 REGRESSION AND CORRELATION COEFFICIENTS FOR GENERATION OF MONTHLY INFLOW VOLUMES 69 7.4.1 DISTRIBUTION OF PROPORTION OF ANNUAL FIRM ENERGY OUTPUT BY MONTH 71 8.1.1 SUBOPTIMIZATION RESULTS, PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 10,000 mkwh . 74 8.1.2 SUBOPTIMIZATION RESULTS, PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 12,000 mkwh . 75 8.1.3 SUBOPTIMIZATION RESULTS, PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 14,000 mkwh . 76 8.1.4 SUBOPTIMIZATION RESULTS, PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 16,000 mkwh 77 8.1.5 SUBOPTIMIZATION RESULTS, PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 18,000 mkwh 78 8.1.6 SUBOPTIMIZATION RESULTS, PORTGAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 20,000 mkwh . 79 vi LIST OF FIGURES FIGURE Page 2.1.1 SEQUENTIAL OR MULTISTAGE SYSTEM. . . . . . . . 10 2.1.2 FULLY INTERACTIVE ARRANGEMENT OF SUBSYSTEMS. . . . . . . . . 10 2.1.3 HIERARCHICAL OR MULTILEVEL DECOMPOSITION OF A SYSTEM . . . . 12 2.2.1 TYPICAL MIXED THERMAL AND HYDROELECTRIC GENERATING SYSTEM. . 14 2.2.2 TWO LEVEL HIERARCHICAL DECOMPOSITION OF A MIXED GENERATING SYSTEM 16 2.2.3 THREE LEVEL HIERARCHICAL DECOMPOSITION OF A MIXED GENERATING SYSTEM 18 6.1 COMPUTER FLOWCHART FOR MAIN PROGRAM. 51 6.2 COMPUTER FLOWCHART FOR DYNAMIC PROGRAMMING SUBROUTINE. . . . 54 6.3 COMPUTER FLOWCHART FOR SIMULATION SUBROUTINE . 57 7.2.1 FREQUENCY DISTRIBUTION FOR ANNUAL INFLOW VOLUME PORTAGE MOUNTAIN SUBSYSTEM . 65 8.1.1 ANNUAL FIRM ENERGY OUTPUT — PRESENT WORTH EXPECTED COST CURVE FOR PORTAGE MOUNTAIN SUBSYSTEM. . . . . . . . . . 80 8.2.1 ANNUAL FIRM ENERGY OUTPUT — PRESENT WORTH EXPECTED COST CURVE FOR PORTAGE MOUNTAIN SUBSYSTEM WITH THERMAL CAPACITY CONSTRAINT 84 8.2.2 CALIBRATED FINITE THERMAL CAPACITY OUTPUT -- COST CURVE. . . 86 8.2.3 RE-CALIBRATED FINITE THERMAL CAPACITY CURVE TO ALLOW FOR DEVIATION FROM INFINITE CAPACITY CURVE . . . 88 9.1.1 ANNUAL FIRM ENERGY OUTPUT — COST CURVES FOR PORTAGE MOUNTAIN AND MICA SUBSYSTEMS 91 9.1.2 MINIMUM COST ALLOCATION OF FIRM ENERGY OUTPUT FOR PORTAGE MOUNTAIN AND MICA SUBSYSTEMS 92 v i i v i i i FIGURE Page 9.1.3 ANNUAL FIRM ENERGY OUTPUT — COST CURVE FOR THE OVERYEAR STORAGE SUBGROUP. 94 9.1.4 MINIMUM COST ANNUAL FIRM ENERGY ALLOCATIONS OVER A RANGE OF OVERYEAR STORAGE SUBGROUP OUTPUTS . „ 95 9.2.1 ANNUAL FIRM ENERGY OUTPUT — COST CURVE FOR SEASONAL AND NON-STORAGE SUBGROUP (Postulated) 97 A C K N O W L E D G E M E N T The author wishes to express a special debt of gratitude to Professor S. 0. Russell for his constant encouragement in this work and the considerable time he expended in many useful discussions. ix CHAPTER 1 INTRODUCTION AND SUMMARY 1.1 GENERAL The energy needs of our society continue to grow and there is no indication that this trend will reverse in the foreseeable future. Since i t is becoming apparent that a l l methods of energy generation have their disadvantages when viewed in the light of their total impact on the environ-ment, i t is vital that we make the best use of existing sources of energy and carefully plan to minimize the adverse effects of new ones. This goal cannot be achieved simply by building the most efficient generating plants; i t is also essential that new plants be tailored to f i t into the existing elec-tric a l system so that the subsequent operation of the whole integrated system is efficient. As the number of generating plants in a system increases and the transmission network develops, a complex interconnected system evolves and the best procedure for efficient integrated operation of the plants becomes increasingly difficult to identify. In a recent (1972) publication [1] Buras states: "The design of multistructure (and multipurpose) systems consisting of several interrelated reservoirs is one of the cardi-nal problems today in water resources engineer-ing. The analysis of such systems with inclu-sion of stochastic elements has hardly been attempted." In this thesis an analytical method is presented for finding an efficient way to operate a hydroelectric generating system which includes 1 2 o n e o r m o r e l a r g e r e s e r v o i r s w h i c h p r o v i d e o v e r y e a r s t o r a g e a n d i s s u b j e c t t o u n c e r t a i n t y r e g a r d i n g i n f l o w s i n t o t h e r e s e r v o i r s i n f u t u r e y e a r s . 1.2 O P E R A T I O N O F H Y D R O E L E C T R I C S Y S T E M S A p r e d o m i n a n t l y s i n g l e p u r p o s e h y d r o e l e c t r i c s y s t e m m a y i n c l u d e r u n o f r i v e r , s e a s o n a l s t o r a g e a n d o v e r y e a r s t o r a g e p r o j e c t s . R u n o f r i v e r p l a n t s h a v e n o s t o r a g e c a p a b i l i t y a n d i n c o m i n g f l o w s m u s t e i t h e r b e u s e d f o r i m m e d i a t e p o w e r g e n e r a t i o n o r b e b y p a s s e d a n d w a s t e d . W i t h s u c h p l a n t s t h e o p e r a t i o n a l p r o c e d u r e i s s t r a i g h t f o r w a r d ; a l l f l o w s a r e s i m p l y u s e d f o r g e n e r a t i o n u p t o t h e l i m i t w h i c h t h e s y s t e m c a n " a b s o r b , a n d t h e r e m a i n d e r a r e b y p a s s e d . S e a s o n a l s t o r a g e p r o j e c t s a l l o w d i s c h a r g e s t o b e d e f e r r e d u n t i l a l a t e r d a t e i n t h e y e a r , b u t , f o r s o m e r e a s o n s u c h a s f l o o d c o n t r o l , a l l t h e w a t e r s t o r e d i s r e -l e a s e d i n t h e c o u r s e o f e a c h y e a r . A l t h o u g h t h e o p e r a t i n g d e c i s i o n c h o i c e i s c o n f i n e d t o t h e r e d i s t r i b u t i o n o f d i s c h a r g e s w i t h i n a p a r t i c u l a r y e a r t h e r e i s u n c e r t a i n t y a b o u t s h o r t - t e r m i n f l o w s , s e a s o n a l t o t a l i n f l o w s , a n d s h o r t - t e r m g e n e r a t i n g l o a d s ; a s a r e s u l t t h e d e t e r m i n a t i o n o f t h e m o s t e f f e c t i v e o p e r a t i o n c a n p r e s e n t a c o m p l e x p r o b l e m . T h e o v e r y e a r s t o r a g e p r o j e c t h a s s u f f i c i e n t s t o r -a g e c a p a c i t y t o p e r m i t w a t e r r e c e i v e d i n h i g h i n f l o w y e a r s t o b e r e t a i n e d a n d u s e d t o g e n e r a t e p o w e r a y e a r o r m o r e l a t e r i f a n d w h e n i n f l o w d e f i c i e n c i e s a r e e x p e r i e n c e d . A l t h o u g h t h e o v e r y e a r s t o r a g e p r o j e c t i s c h a r a c t e r i s t i c a l l y l a r g e a n d l e s s s e n s i t i v e t o s h o r t - t e r m f l u c t u a t i o n s i n i n f l o w s a n d g e n e r a t i n g l o a d s t h a n t h e s e a s o n a l s t o r a g e p r o j e c t , t h e r e i s i n t e r d e p e n d e n c y b e t w e e n o p e r a t i o n s f r o m o n e y e a r t o a n o t h e r b u t n o p o s s i b i l i t y o f f o r e c a s t s o f f u t u r e y e a r ' s i n -f l o w s . E f f i c i e n t o p e r a t i o n o f a n a l l h y d r o s y s t e m w h i c h i n c l u d e s a n u m b e r o f r u n o f r i v e r p l a n t s a n d a s i n g l e s e a s o n a l s t o r a g e p r o j e c t i s a c h i e v e d w h e n t h e 3 run of river plants generate at their maximum capability (determined by the natural river flows) and the storage plant only generates to meet the residual load. In the case where a system consists of a number of run of river plants and thermal plants the run of river plants would again be operated at their maximum capability and then the thermal plants would be called upon, the most efficient f i r s t and then progressively the less efficient, until the residual load is met. In this way the overall efficiency of the generating system would be maximized as the total load would be met with the minimum con-sumption of fuel by the thermal plants. For a system which includes thermal generating plants, possibly run of river plants, and one or more hydroelectric plants with seasonal storage than the most efficient operation cannot be achieved by a simple policy which fav-ours generation from one type of plant before calling upon the others. There are now many possible ways in which the residual load, over and above that which can be met by the run of river plants, can be met depending on the timing of the release of the stored water. If the natural inflows can be predicted with reasonable accuracy, as is the case for many snowmelt fed storage hydro projects in Canada, then the operating problem is one of sequential decision making under deterministic conditions and can be readily analyzed by determinis-tic dynamic programming. The analysis is confined to a period of one year or less but this is sufficient to produce optimal results since operation of this system in one year does not influence operation in subsequent years. The Aluminum Company of Canada's generating system in Quebec [2] and the British Columbia Hydro and Power Authority's Campbell River system are two examples of seasonal storage hydro systems whose operating decisions are being determined by this dynamic programming approach. 4 The existence of overyear storage in a hydroelectric system intro-duces additional difficulties to the operating decision analysis. Since a wide range of storage volumes are now possible at the end of any year the most desireable year end conditions (such as empty) cannot be specified in advance as was the case for the seasonal storage project. The best operating decisions during the course of a particular year cannot therefore be determined until one can determine how the volume left in storage at the end of the year will affect operation in the years following. This influence of one year's operation on another introduces some additional complications which do not arise when determining the optimal short-term operating policies for generating systems without overyear storage capabilities as discussed above. The interdependency of operation in one year with operation in follow-ing years extends indefinitely into the future; hence the decision process has an unbounded or indefinite time horizon. Furthermore, although inflow fore-casting is generally feasible during a current runoff year, i t is not possible to forecast inflows beyond the current year and therefore uncertainty regarding the outcomes of storage decisions which are inflow dependent cannot be avoided. The difficulties of analysing the operation of a large generating system due to its inherent complexity, both in conjunction with the formulation and the compu-tation, are seriously compounded by these additional complications of unbounded time horizon and uncertainty. A variety of techniques for resolving these difficulties have been demonstrated in Water Resources applications. Uncertainty and the unbounded time horizon may be avoided i f the analysis is confined to a finite time period and uncertainty eliminated by considering a deterministic hydrology representing either a c r i t i c a l inflow sequence [3] or, alternatively, by repeating the 5 analysis a number of times over a set of simulated typical inflow patterns [4]. The optimal policies obtained from this form of analysis provide valu-able insight into the operation of the system, however, they relate to an a r t i f i c i a l deterministic situation and not to the reality of the situation in actual operation where operating commitments must be made in the face of uncertainty. A more direct probabilistic approach to deal with uncertainty is the use of Markov models which, provided certain statistical requirements are met, provides a computationally attractive approach. A variety of Markov models used in the analysis of the regulation of a single reservoir have been re-viewed by Gablinger and Loucks [5] and Buras [1]. The problem of overall complexity when analysing a large generating system over a substantial time period has been resolved by decomposing the overall problem into a set of sub-problems, which individually are more amen-able to analysis, and then performing some form of reco-ordination based on the results obtained. Hall, Butcher, Esogbue [6] [7] have demonstrated an optimisation procedure of this type based upon Linear Dynamic Decomposition Programming [8] for a multireservoir multipurpose water resource system but using a deterministic inflow sequence. Buras [1], Hall [9], and Chow [10] have described more general forms of decomposition approaches where separate optimisation procedures are used in the analysis of the subsystems and in the reco-ordination of these subsystem results. Although this approach does not ensure that a true global optimum is obtained for the overall system i t permits investigation of considerably more complex systems than can be con-sidered by rigorous mathematical methods and introduces a high degree of organization to the system components. 6 The problem considered in this study, i.e., the long-term operation of a large hydrothermal generating system with overyear storage project(s), embodies these difficulties of unbounded time horizon, uncertainty, and system complexity. The procedure developed for the analysis of this operation was designed to resolve a l l three of these difficulties. The system was first decomposed into a set of interconnected subsystems having a hierarchical control structure which retains the essential complementary and interractive aspects of the generating units in the context of long term firm energy production. A Markov model was adopted to represent the long term annual operation of an individual subsystem comprising of an overyear storage hydro project and a thermal plant. Howards method for the analysis of Markov decision processes was used in the determination of the optimal operating policy for this subsystem. A departure from Howard's basic algorithm was introduced by using a 12-stage deterministic Dynamic Programming analysis to optimize the within-year operation of the subsystem and simultaneously complete the Policy Improvement phase of the algorithm. The results of this subsystem optimization, i.e., suboptimization, provide a measure of the cost of long-term operation as well as monthly generating discharges for the hydro project. An example of the application of this method is given for the British Columbia Hydro system. This system consists of a number of small to medium-sized hydroelectric plants of limited or zero storage, a 900MW thermal plant at Burrard, and a 2200MW overyear storage hydro project with a live storage capacity of 15 million acre f t . at Portage Mountain on the Peace River. A second overyear storage hydro project at Mica Creek on the Columbia River is currently nearing completion and will have a generating capacity of 2600MW and live storage capacity of 12 million acre f t . 7 The long-term operation of the subsystem, comprising of the Portage Mountain project and an allocation of thermal plant capacity, was analyzed to determine the optimal (i.e., minimum thermal energy consumption) discharge policy for a prescribed annual firm energy output. The suboptimization pro-cedure was found to be computationally efficient and was repeated for this sub-system over a range of thermal capacity allocations and subsystem outputs. These results were used to produce an output/cost curve which reveals the poten-t i a l benefits which may be obtained in integrated operation of the hydro and thermal plants represented within the subsystem. This curve was then used as the basis of a procedure for co-ordinating the long-term operation of the Portage Mountain and Mica Creek projects. This co-ordinating procedure was finally extended to the co-ordination of these two projects with the remainder of the generating system. CHAPTER 2 2.1 OPTIMIZATION OF COMPLEX SYSTEMS Almost any otpimization problem can be solved, at least in principle, by existing mathematical methods but computational difficulties arise when the dimensionality of the problem becomes large; this is generally the case when large or complex systems are considered. These computational difficulties are overcome to some extent when the overall optimization problem is reduced to a set of subproblems which are dimensionally small enough to permit the optimi-zation of each independently of the remaining subproblems. In this manner com-plex systems can be more easily handled. Methods which reduce the dimensionality by setting up subproblems each of lower dimensionality are called decomposition techniques. These methods do not prescribe the optimization procedure to be used in the subproblems but pro-vide a framework for subdividing or partitioning the original problem and sub-sequently recombining the solved subproclems. No general decomposition approach nas yet been developed for dealing with complex systems but due to the practical importance of optimizing such systems i t is the subject of active investigation. The partitioning of a system may be accomplished in many ways but the combined results of optimizing the subproblems individually will not reflect the optimal behaviour of the overall system unless interraction of the subsystems is taken into account. The decomposition approach adopted in solving a particular system problem must therefore not only reduce the total system problem to a set of subproblems of manageable size but also provide a feasible basis for 8 9 c o - o r d i n a t i n g t h e s e s u b p r o b l e m s . T h e s e t w o r e q u i r e m e n t s a r e , o f c o u r s e , i n t e r -r e l a t e d a n d o f t e n c o n f l i c t i n g . R e d u c t i o n i n t h e s i z e o f t h e s u b p r o b l e m s s i m p l i -f i e s t h e m b u t i n c r e a s e s t h e n u m b e r o f s u b p r o b l e m s n e e d e d t o r e p r e s e n t t h e o v e r -a l l s y s t e m t h u s i n c r e a s i n g t h e c o m p l e x i t y o f t h e c o - o r d i n a t i o n p r o b l e m . I n t h e f i e j d o f s y s t e m s e n g i n e e r i n g , w h i c h d e a l s w i t h t h e d e s i g n a n d c o n t r o l o f p r o p o s e d o r e x i s t i n g s y s t e m s , a f u r t h e r r e q u i r e m e n t a r i s e s . I f t h e t h e o r e t i c a l s y s t e m p e r f o r m a n c e i s t o b e a c h i e v e d i n p r a c t i c e , i t i s e s s e n -t i a l t h a t t h e r e s u l t s o f t h e d e c o m p o s i t i o n - r e c o - o r d i n a t i o n p r o c e s s o f a n a l y s i s b e s u i t a b l e f o r i m p l e m e n t a t i o n i n t h e f o r m o f d e s i g n o r c o n t r o l d e c i s i o n s a p p l i c a b l e t o t h e r e a l w o r l d s y s t e m . W h e n a s y s t e m i s d e s c r i b e d i n t e r m s o f a s e t o f s u b p r o b l e m s o r s u b -s y s t e m s w i t h i n t e r r a c t i v e r e l a t i o n s h i p s b e t w e e n t h e s u b s y s t e m s r e p r e s e n t e d b y d i r e c t e d l i n k s , t h e n t h e r e s u l t i n g c o m p l e x r e v e a l s t h e s t r u c t u r e o f t h e s y s t e m i n i t s d e c o m p o s e d f o r m . I f t h e s u b s y s t e m s a r e , f o r e x a m p l e , l i n k e d i n a s e r i e s a s s h o w n i n F i g u r e 2.1.1, w h e r e d e c i s i o n s a p p l i e d i n o n e s u b s y s t e m o n l y a f f e c t s u b s y s t e m s t o t h e r i g h t , t h e n t h e o v e r a l l s y s t e m h a s a s e q u e n t i a l s t r u c t u r e . P r o v i d e d t h a t r e q u i r e m e n t s o f s e p a r a b i l i t y [11] o f t h e r e t u r n s f r o m e a c h s u b -s y s t e m i n t h e s e q u e n c e a r e m e t , a s y s t e m w h i c h c a n b e d e c o m p o s e d i n t h i s f a s h i o n c a n b e o p t i m i z e d u s i n g a s t a n d a r d D y n a m i c P r o g r a m m i n g t e c h n i q u e . M o r e g e n e r a l l y , h o w e v e r , t h e s u b s y s t e m s i n t e r a c t s i m u l t a n e o u s l y a n d e a c h i n d i v i d u a l s u b s y s t e m i n t e r a c t s w i t h a n u m b e r o f o t h e r s u b s y s t e m s . I n s u c h s i t u a t i o n s a d i a g r a m o f t h e s u b s y s t e m s a n d t h e i r l i n k a g e s m a y n o t r e v e a l a n y p a r t i c u l a r s t r u c t u r e w h i c h c o u l d b e e x p l o i t e d t o f a c i l i t a t e a n a l y s i s . A d e -c o m p o s e d s y s t e m o f t h i s t y p e i s s h o w n i n F i g u r e 2.1.2. A n a l y s i s o f a s y s t e m h a v i n g t h i s t y p e o f s t r u c t u r e p r e s e n t s p r o b l e m s o f a c o m b i n a t o r i a l n a t u r e a n d a p r o h i b i t i v e c o m p u t a t i o n a l l o a d . U n d e r t h e s e c i r c u m s t a n c e s a s i m p l i f i e d f o r m 10 S u b s y s t A S u b s y s t B > 1 F I G U R E 2.1.1 S E Q U E N T I A L O R M U L T I S T A G E S Y S T E M F I G U R E 2.1.2 F U L L Y I N T E R A C T I V E A R R A N G E M E N T O F S U B S Y S T E M S 11 structure must be imposed upon the decomposed system to permit an analysis to be made. The adoption of a hierarchical or multilevel structure has been suggested in many fields of systems investigation to facilitate the analysis without losing the essential interractive character of the system components [12]. The hierarchical structure breaks the decision making or control pro-cess into a number of levels. (See Figure 2.1.3). At the lowest level decis-ions are made for an individual subsystem independently of the other subsystems but subject to some form of operating policy imposed by the next higher level in the hierarchy. This policy ensures that a subsystem operates in a co-ordi-nated fashion with respect to other subsystems which are members of the same subgroup; the subgroup containing a set of subsystems which are intensely inter-active. At progressively higher levels in the hierarchy the subgroups, and thus subgroup policies, embrace increasing numbers of the subsystems until policy made at the highest level influences, via intermediate levels of the hierarchy, a l l of the subsystems. No direct interconnection occurs between subsystems or subgroups on the same level so that, although the hierarchical structure resembles a sequential decision process with branching, decisions made at one level of the hierarchy may affect both higher and lower levels in the sequence(s). (See Figure 2.1.3). It is unlike the sequential decision process shown in Figure 2.1.1 where a subsystem could only influence subsys-tems following in the sequence but not preceding subsystems. The hierarchical approach has been suggested in connection with the optimization of short-term power dispatching in a generating system [13]. Al-though the nature of this particular problem differs considerably from optimi-zation of the long-term operation of a hydro thermal generating system considered S y s t e m C o o r d i n a t i n g L e v e l Subsystem Operat ing L e v e l S u b g r o u p C o o r d i n a t i n g L e v e l Subsystem] A It D E FIGURE 2.1.3 HIERARCHICAL OR MULTILEVEL DECOMPOSITION OF A SYSTEM 13 in this thesis, the hierarchical approach is well suited to many types of control problems and was adopted in this study. Following the adoption of the hierarchical approach to the analysis of a system involving a number of interacting components there are many possibilities for the decomposition to subsystems and the grouping of sub-systems at each level of the hierarchy. These possibilities must be assessed in terms of meeting the requirements discussed above. These requirements may be summarized as: (a) the complexity of the individual subsystem analysis must be within the capability of existing decision analysis techniques as well as computational capabilities; (b) vital interractive aspects of the subsystem operation must not be unduly suppressed or lost; (c) co-ordination of the subsystems and subgroups must be feasible; (d) the results must be produced in a form which is implementable in subsequent operation of the real system. The implication of these requirements are considered in the following section in connection with some hierarchical decomposition possibilities for the analysis of the long-term operation of a mixed hydro and thermal generating system with overyear storage. 2.2 DECOMPOSITION OF A MIXED HYDRO AND THERMAL GENERATING SYSTEM IN THE FORM OF A HIERARCHICAL SYSTEM The elements of a typical hydro and thermal generating system are shown schematically in Figure 2.2.1. The system consists of hydro electric projects having overyear storage, seasonal storage, and run of river plants with no storage as well as a number of thermal plants. A network of transmis-sion lines connects the generating units to the load centre and can be considered w S e a s o n a l S to rage Hydro Project O v e r y e a r S to rage Hydro Pro jec t Run of R i v e r ^Hydro Pro jec t Ruif 'of /^p7\ R iver tf Hydro P r o j e c t T h e r m a l P l a n t Run v o f R i ve r ^ Hydro Project [ mil O v e r y e a r S t o r a g e Hydro P r o j e c t FIGURE 2.2.1 TYPICAL MIXED THERMAL AND HYDROELECTRIC GENERATING SYSTEM 15 to provide f u l l interconnection between a l l components of the system. The transmission network arrangement and capacity will be determined by the need to meet peak loads and other short-term requirements. Long-term operation based on monthly and annual energy flows will not normally be a determining factor in the design of the transmission network and conversely the trans-mission system will not be a constraining factor in establishing the long-term operating policy for generation. Full interconnection of the generating components results in a sys-tem which, from the point of view of its operation to meet a system load, is fully interactive. The interaction occurs simultaneously in both direc-tions along each link. The optimum energy contribution of one generating unit over a period of time will depend upon operating conditions experienced by that generating unit during the period as well as the operating condition of a l l of the remaining generating units in the system. An additional compli-cation is introduced by the presence of storage in the system; operation in one time period will affect operating conditions in future time periods. The fully interconnected structure of the generating system resembles that shown in Figure 2.1.2 and discussed in Section 2.1. Direct analysis of the long-term operation of such a system is not possible. Analysis and optimi-zation of the total system operation can, however, be achieved by reorganizing the system generating components into a set of subsystems, selecting suitable subgroups of subsystems, and then imposing a hierarchical structure on the system. Figure 2.2.2 demonstrates one conceivable decomposition of the system but one which fails to meet the requirements discussed in Section 2.1. Here a subsystem contains a l l generating units with similar operating characteristics, 16 A Coordination Level for Subsystems Storage Hydro Projects Thermal Subsystem Non Storage Subsystem Hydro Subsystem FIGURE 2.2.2 TWO LEVEL HIERARCHICAL DECOMPOSITION OF A MIXED GENERATING SYSTEM 17 e.g., a l l hydro projects with storage are assigned to one subsystem. Opti-mized operation of the individual subsystems in this case does not reflect interaction between generating units which have very different, and there-fore potentially compementary, operating characteristics. Interactions of this type are only considered at higher levels in the decision hierarchy after considerable aggregation of individual plant characteristics has taken place in the analysis of each subsystem. In addition, the simultaneous optimi-zation of a group of storage hydroelectric projects within a subsystem raises problems of dimensionality and hence computational difficulties. This is par-ticularly true when the random properties of the individual inflows are taken into consideration. Arventidas [14] has demonstrated a method to overcome the dimensionality problem by using potential energy as a composite measure of the storage in a number of reservoirs. The use of a composite state variable pre-supposes the amount of energy extractable per unit of storage. Although i t provides useful results in the form of optimum composite policies this method does not produce operating decisions which are appropriate for implementation at the individual storage project level. A more suitable decomposition of the generating system is shown in Figure 2.2.3, and was adopted in this study. It provides for f u l l interaction between storage hydro and thermal plants by introducing a hypothetical thermal plant within each subsystem which contains a storage hydro plant. Since the system is fully interconnected the location of the thermal plants in the de-composed system, or a subdivision of them, is immaterial provided that the total thermal capacity distributed to the various subsystems is equal to the actual total thermal capacity available in the system. Hypothetical thermal plants are not included in subsystems containing non-storage hydro plants as these ^ / \ ^ System Coordinating Level Overyear Storage Overyear Storage Hydrothermal Hydrothermal Subsystem I Subsystem 1 Residual Thermal Subsystem Seasonal Storage Hydrothermal Subsystems Run of River Subsystems ElGURE 2.2.3 THREE LEVEL HIERARCHICAL DECOMPOSITION OF A MIXED GENERATING SYSTEM 03 19 do not offer any significant long-term operational benefits as a result of integration. Remaining thermal plant capacity which is not assigned to integrated operation is included in a single purely thermal subsystem. At the second level of the hierarchy the subgroups include sub-systems which interact most strongly with respect to their long-term operating policies to the system firm energy output. One subgroup includes a l l of the overyear storage subsystems but excludes seasonal storage, run of river, and residual thermal subsystems. In the case when overyear storage hydro projects provide a major contribution to the system firm energy out-put, the role of projects where operation in one year cannot be made to influence their operation in subsequent years is confined to modifying the shape of the annual firm energy load imposed on the overyear storage sub-systems (see Chapter 7). The overyear storage subgroup was the principal concern of this study. The decomposition of the remaining subgroups and their co-ordination in the context of overall long-term system policy is discussed in Chapter 9. 2.3 SUBOPTIMIZATION OF OVERYEAR STORAGE SUBSYSTEMS An analysis to determine the long-term optimal operating policy for an overyear storage subsystem can be performed once a firm energy output and thermal capacity allocation are prescribed. As this information is not avail-able prior to completing the total system analysis i t is necessary to treat sub-system energy output and thermal capacity allocation parametrically and perform the analysis over a range of values for each of these parameters [9], [1]. The results of such an analysis provides a relationship between a policy imposed on a storage subsystem from the second level in the hierarchy (in terms of an energy 20 demand and a thermal capacity allocation) and the subsystem output response (in terms of the cost of meeting the policy measured in thermal energy con-sumption) • The character of the decision process in the long-term operation of an overyear storage subsystem and its implication in the suboptimization is discussed in Chapter 3. 2.4 CO-ORDINATION OF OVERYEAR STORAGE SUBSYSTEMS The objective to be achieved in the co-ordination of the subgroup of overyear storage subsystems is to extract the maximum firm energy contri-bution of the subgroup to the overall system firm energy demand. This ob-jective is to be met subject to two principal constraints. The incremental or marginal cost of producing firm energy in this subgroup must be no greater than the cost of generating additional firm energy elsewhere in the system, and the firm energy production of the entire generating system must meet the demand. The selection of a suitable co-ordination procedure is dependent upon the characteristics of the subsystems revealed by suboptimization. The co-ordination procedure adopted for the overyear storage subgroups is described in Chapter 9. C H A P T E R 3 L O N G - T E R M O P E R A T I O N O F T H E H Y D R O T H E R M A L S U B S Y S T E M 3.1 I N T E G R A T E D O P E R A T I O N O F T H E R M A L A N D H Y D R O G E N E R A T I N G P L A N T S H y d r o e l e c t r i c a n d t h e r m a l g e n e r a t i n g p l a n t s h a v e c o n t r a s t i n g e c o n o m i c a n d t e c h n i c a l c h a r a c t e r i s t i c s . T h e c o m b i n e d o p e r a t i o n o f t h e t w o t y p e s o f p l a n t w i t h i n t h e h y d r o t h e r m a l s u b s y s t e m p r o v i d e s a n o p p o r t u n i t y t o e x p l o i t t h e i r c o m p l e m e n t a r y c h a r a c t e r i s t i c s a n d a c h i e v e t h e m o s t f a v o u r a b l e o v e r a l l p e r f o r m a n c e . A h y d r o e l e c t r i c p r o j e c t i n v o l v e s a l a r g e c a p i t a l c o s t c o m m i t m e n t ( i . e . , f i x e d c o s t ) i n i t s i n i t i a l c o n s t r u c t i o n b u t l o w s u b s e q u e n t o p e r a t i n g c o s t s ( i . e . , v a r i a b l e c o s t ) w h i c h a r e v i r t u a l l y i n d e p e n d e n t o f t h e a m o u n t o f e l e c t r i c a l e n e r g y p r o d u c e d . A t h e r m a l p l a n t i n v o l v e s , i n c o m p a r i s o n , a p r o p o r t i o n a l l y l o w c a p i t a l o r f i x e d c o s t b u t h i g h o p e r a t i n g c o s t s . T h e s e r u n n i n g c o s t s a r e r e l a t e d d i r e c t l y t o t h e f u e l c o n s u m e d a n d t h e r e f o r e t o t h e a m o u n t o f e l e c t r i c a l e n e r g y p r o d u c e d . T h u s w h i l e a h y d r o p l a n t o f f e r s v i r t u a l l y n o o p p o r t u n i t y f o r a n y m a j o r e c o n o m i c d e c i s i o n m a k i n g f o l l o w i n g i t s i n i t i a l c o n s t r u c t i o n t h e r e i s a c o n t i n u i n g o p p o r t u n i t y t o c o n t r o l t h e e c o n o m i c c o s t s a s s o c i a t e d w i t h t h e t h e r m a l p l a n t . F o r a h y d r o t h e r m a l s u b s y s t e m c o n s i s t i n g o n l y o f a n e x i s t i n g o r c o m m i t t e d h y d r o p r o j e c t a n d a n a l l o c a t i o n o f e x i s t i n g t h e r m a l c a p a c i t y t h e n t h e o n l y s i g n i f i c a n t f o r m o f e c o n o m i c c o n t r o l o r d e c i s i o n m a k i n g a s s o c i a t e d w i t h i t s s u b s e q u e n t o p e r a t i o n i s i n t h e r m a l f u e l c o n s u m p t i o n a n d h e n c e t h e r m a l g e n e r a t i o n . T h e h y d r o p r o j e c t r e l i e s u p o n a n u n c o n t r o l l e d e n e r g y s u p p l y i n t h e f o r m o f t h e n a t u r a l i n f l o w t o t h e r e s e r v o i r . T h i s i n f l o w i s o f t e n s u b j e c t t o 21 22 large variations from year to year. The reservoir is capable of storing this incoming energy subject to the limitations of its storage capacity. The energy extracted per unit volume of discharge through the turbines is not fixed but varies directly with the turbine head and thus the stored energy at the time of its discharge. This uncontrolled and variable energy source and storage dependent energy conversion process contrasts with the controllable energy supply (assuming adequate thermal fuel sources are available) and un-varying energy conversion process of the thermal plant. In the combined operation of thermal and hydroelectric plants the dependability of the thermal energy source may be used to make up deficiencies in the output of the hydro plant resulting from deficient inflows and/or storage depletion. This alone, however, is a myopic operating policy which overlooks the greater benefits which may be obtained by integration over a longer period of time. The thermal plant can be used at non-critical times to relieve the hydro plant of load and hence, by retaining the inflow in storage, raise the reservoir level. As a result the hydro plant can subsequently operate at higher heads thereby extracting more electrical energy per unit volume of water discharged. In effect this procedure commits economic costs at one point in time, in the form of increased thermal fuel consumption, in'order to achieve future benefits in the form of enhanced hydro generation and therefore u l t i -mately net savings in thermal fuel costs. The benefits resulting from load relief operations can be substantially greater than the cost of thermal fuel involved. While the cost of load relief is proportional to the increase in storage volume, the resulting increase in energy extraction per unit discharge may be sustained over a considerably larger dis-charge volume. (An extreme example of this aspect is in the i n i t i a l f i l l i n g of 2 3 a hydro reservoir where the i n i t i a l dead storage sustains operating heads for the entire l i f e of the project). Furthermore, as the energy invested in load relief i s in effect stored in the form of increased reservoir stor-age then a large proportion of this energy is recoverable (although after a considerably longer period of time than is usually associated with pumped storage cycles). The head increase attributable to a load relief decision may be reduced or eliminated by drawing down the reservoir 'or by the arrival of inflows which exceed the available storage capacity and result in spilling. Thus load relief benefits are dependent upon the inflows which occur follow-ing the load relief decision and also upon subsequent operating decisions. The characteristic length of time over which reservoir inflows will influence the benefits obtained from a load relief decision will depend upon the storage capability of the reservoir. A hydro project with seasonal stor-age would not be capable of sustaining boosted reservoir levels beyond the year in which the load relief decision was made. In the case of the overyear storage project i t is possible that the increased reservoir level obtained from load relief decisions in the overyear storage subsystem can only be ..assessed by considering the consequences of such a decision over a number of years of operation, and, as has already been discussed, the consequences will be de-pendent upon the inflows occurring and subsequent operating decisions during this period. Furthermore, as operation in any particular year will affect oper-ation in a following year which in turn affects operation in subsequent years, then the operational decision process cannot be logically terminated at- any point in time and must be regarded as a decision process having an indefinite time horizon. 24 The implications of the above factors and their resolution in the suboptimization process are discussed in the remainder ofwthis chapter. Section 3.3 discusses the predictability of inflows to the overyear stor-age reservoir, Sections 3.3. and 3.4 describe how the problems of indefinite time horizon and uncertainty of future inflows were resolved. 3.2 INFLOW PREDICTION FOR THE OVERYEAR STORAGE HYDRO PROJECT The large storage reservoirs considered in this study are fed pri-marily from snowmelt and i t is usually possible to provide a relatively accur-ate forecast of the total annual inflow volume during a current runoff year. The monthly inflow pattern over the current year is also predictable with the major portion of the total annual inflow occurring during the early summer months. The combination of an accurate forecast and the insensitivity of very large reservoirs to minor forecast errors [22] allows for virtually determinis-tic operation of the reservoir during this period (provided a year end goal can be specified). It is not, however, possible to make prediction with meaningful accur-acy as to what total inflow volume will occur in a year following the current runoff year. The historic record of annual inflows provides information on its probabilistic nature but provides no evidence of any significant serial corre-lation from one year to another (see Chapter 7). The future annual inflow volumes were therefore considered as independent random quantities. As i t is necessary, with the existence of overyear storage, to assess load relief decisions beyond the current runoff year and this assessment is part-ly dependent upon inflows both in the current and future years, then the assess-ment will be made under certainty during the current year and under uncertainty in the years following. 25 3 . 3 THE INDEFINITE TIME HORIZON The operation of a hydrothermal subsystem which includes an over-year storage hydro project is a sequential decision process which continues indefinitely into the future and which has a considerable degree of interdepend-ency between years. It is unlikely that the analysis of any long-term or infinite operational process would be used to provide firm operational commitments beyond the f i r s t few years. One would expect that operating decision commitments (as opposed to contractual commitments) beyond this time would only be made follow-ing a reappraisal of the situation. It i s , however, necessary to consider the operational decision process beyond the commitment term to avoid arriving at short sighted policies which may place later operation in jeopardy. The optimization of a process with an indefinite or unbounded time horizon is not possible by simply maximizing the sum of the sequential returns as this would involve comparison of infinite sums of returns. This problem can be overcome in three ways: by optimizing the average return over some fixed time period, e.g., one year; by discounting the returns and optimizing their combined present worth; or by terminating the analysis at some time in the future and assigning suitable values or penalties depending on the conditions at the terminal point in time [15]. The choice of method for resolving the indefinite time horizon depends upon the nature of the problem being analyzed. The average gain per year criterion will reflect th e infinite nature of the process and as a result is insensitive to large i f infrequent losses. This can be a severe drawback whenever sustained poor performance over a few time periods is unacceptable or can lead to the termination of the enterprise (e.g., gamblers ruin). When only the i n i t i a l decisions in the un-bounded decision sequence are in fact firm commitments an optimal policy based 2 6 on average gain per year may result in negligible returns or even losses over the i n i t i a l stages yet from the decision maker's point of view, these may well be the most important stages for performance evaluation . Generally the average gain per stage criterion is most suited to processes involving decision commitments over a large number of stage repetitions and was con-sidered inappropriate in the analysis of annual operating policy for the hydrothermal subsystem where good performance is essential over both the immediate short-term and long-term future. The terminal approach is a computationally attractive method of resolving the indefinite horizon problem. It requires a value for each of the set of possible discrete terminal states (e.g., in the case of the over-year storage subsystem the volume of water stored) which reflects the value of each state to future optimal operation. These terminal state values are used to ensure that optimal operation prior to the terminal stage is compatible with optimal operation following i t . Although this introduces the useful con-cept of decomposition of the infinite stage process into two parts, before the terminal stage and after i t , and thereby permits a finite stage analysis to be performed up to the terminal stage, i t does not fully resolve the difficulty. It is s t i l l necessary to obtain the terminal state values and these are s t i l l dependent upon the indefinite time horizon process following the terminal stage. In determining the present worth of a stream of returns over time the discount factor progressively (i.e., geometrically) unweights the returns as they occur further into the future. The unweighting operation is appropri-ate in a non-monetary sense as well as the more usual financial one since i t reflects the diminishing commitment and increasing doubt associated with decisions 27 and returns which are further away in time. For problems with an unbounded time horizon the present worth method provides a rational basis for the analy-sis and lies between arbitrary analysis of short-term operation approach of considering short-term operation only and the average gain per stage approach to the infinite stage process; i t escapes the myopia of the former and the insensitivity to short-term failure of the latter. Its position in the spec-trum between these two extremes depends on the discount factor chosen. Financial interest rates, and hence discount factors, provide a guide to the magnitude of discount factors both in monetary and non-monetary situations. Non-monetary theories of interest [16] suggest that interest rates under minimal risk conditions are determined by human time preferences when time and money interact and represent an aggregate measure of human time per-spective. Although a precise value of the discount factor cannot be prescribed for the analysis of an indefinite time horizon decision process experience has shown that the optimal policies obtained are relatively insensitive to minor variations in the discount factor over the range of typical monetary values. In this study an annual discount factor of 0.926 (equivalent to an interest rate of eight per cent) was adopted for most of the computations. Sensitivity analysis showed that in the range of annual discount factors from 0.952 to 0.894 there was l i t t l e change in the optimal policies. 3.4 ASSESSMENT OF LOAD RELIEF DECISIONS UNDER UNCERTAINTY When there is uncertainty as to the inflow volumes and hence out-come following a load relief decision but the probability of each particular outcome occurring and its associated returns are known, then the mathematical expectation of return provides a single valued measure of the return for the decision. Although this is a computationally convenient measure of return under 28 uncertainty i t suffers two major drawbacks. Expectation, by definition, provides a measure of the average out-come following many applications of that decision under identical circumstan-ces. It may never reflect the actual return following a single implementation of that decision. Where the decision is enacted infrequently, e.g., annually, then a large number of repetitions would involve an inordinate period of time before the average return for a decision, given the same circumstances at each repetition, was approached. Generally a decision maker will be more concerned about the return from, say, an annual decision over the course of a few years rather than over the course of many decades. A second drawback involves the possibility of cri t i c a l or disastrous outcomes which result in extremely high losses and perhaps terminate the oper-ation. Situations such as these would often be associated with very low prob-abilities. In this case the product of the large cost but low probability may not reflect the true implication to future operation i f , for example, further operation is not possible. This shortcoming was avoided in this study by assigning such large arbitrary losses to these c r i t i c a l outcomes so that even i f a decision causes one to occur with the lowest discrete prob-ability adopted in the analysis, then the expected cost will be high enough to ensure that the decision will be rejected. The use of a u t i l i t y curve of this type and resulting externalization of c r i t i c a l or disastrous events (such as a number of consecutive years of drought) from the analysis can be justified on the grounds that unconventional or extreme measures would be taken in the face of such events and these measures are outside the scope of an analysis of normal operation involving only conventional decisions. 29 3.5 SUBOPTIMIZATION SCHEME FOR THE HYDROTHERMAL SUBSYSTEM The formulation for suboptimization of the integrated operation of the hydrothermal subsystem at a prescribed firm energy output level is broken down into two parts: the analysis of subsystem operation during any current year and analysis of the long-term operation over subsequent years. In the analysis of the current year's operation the total annual in-flow volume and its monthly distribution is treated deterministically. The problem is formulated in Chapter 4 using discrete Dynamic Programming with 12 monthly stages. This determines the annual cost of operating the subsys-tem and the year end reservoir level for a given i n i t i a l reservoir level and a given set of year end reservoir state values. The long-term operation of the subsystem is treated as a Markov decision process and the optimization objective is the minimization of the present worth of the future expected costs. This is formulated in Chapter 5. This phase of the analysis utilizes the annual operating costs and reservoir level changes determined from the analysis of the annual operation of the sub-system. CHAPTER 4 ANALYSIS OF ONE YEAR OPERATION OF THE SUBSYSTEM In this chapter the optimal operation of the hydrothermal subsystem over a "current" runoff year is determined using Dynamic Programming. The form-ulation and solution of similar problems using Dynamic Programming has been demonstrated in a number of publications [9], t l ] . The subsystem under consideration consists of an overyear storage hydro project with known storage capabilities and known generating characteristics. It also includes an allocation of thermal plant generating capacity. The oper-ational objective for the year is to meet a specified total annual firm energy output, prescribed in the form of 12 monthly energy demands, with a given total annual inflow. This objective is to be achieved within the long-term objective of consuming the minimum amount of thermal energy which represents the only signi-ficant variable operating cost in the subsystem (see Section 3.1). A number of practical constraints are present in the operation of the subsystem, e.g., reservoir level limits, generating plant limits, etc., and are readily incorporated into the Dynamic Programming analysis. For purposes of clar-ity these will be excluded from the formulation in this chapter and in Chapter 5. The detailed operating constraints used in the computation are described in Chap-ter 6. 4.1 DYNAMIC PROGRAMMING FORMULATION Operating during the year is broken down into 12 monthly intervals which constitute the stages of the Dynamic Programming analysis. The reservoir elevation at the beginning of a one-month stage represents the (discrete) state 30 31 variable. The following notation is used: m = month, m = 1,. . .,12; i = a discrete reservoir elevation (state) at beginning of month m, i =1/ 2,. . .N where N is the number o? discrete states; S(i ) = reservoir storage volume when reservoir elevation is i ; m d = a hydro plant generating discharge volume decision during month m; I = the reservoir inflow volume in month m; m F = the firm energy demand imposed on the sub-system during month m. The continuity equations S(i ,,) = S(i ) + I - d 4.1 m+1 m m m is used to determine the storage volume at the beginning of month m+1 when the elevation at the beginning of month m was im and during the month m there occurred an inflow 1^ and a volume d^ was discharged for generation purposes. Equation (4.1) is used in conjunction with the reservoir storage/elevation curve to determine the reservoir elevation at the beginning of stage m+1. This is represented by the state transformation equation: i + 1 = t ( i , d , I ) 4.2 m i m m m The electrical energy generated by the hydroelectric plant will be dependent upon the discharge volume and the turbine head which in turn depends upon the reservoir elevation. The reservoir elevation will not gen-erally remain constant over the monthly stage, however in the case of large overyear storage reservoirs changes in reservoir elevation occur relatively 32 smoothly. The arithmetic mean of the reservoir elevations at the begin-ning and end of a one-month stage provides an acceptably accurate esti-mate of the mean monthly head for the purposes of long-term firm energy calculations. Thus the hydro energy generated during the course of a month will be dependent upon d , i and i ,,. As i ,, is determined c c m m m+1 m+1 from Equation C4.2) and is in turn dependent upon i , d and I , then the m m m energy conversion function W(i , d , I } is defined as a function of m m m these three variables. If the generation in month m is E then: m E = W(i , d , I ) 4.3 m m m m The cost of meeting the subsystem firm energy demand for the month m is a function of the deficiency between hydro output for the month and the firm energy demand. Thus: C = R(F - E ) 1 4.4 m m m In order to avoid hydro generation in excess of the firm energy demand, and hence negative costs, the constraint W(i , d , I ) < F m m m — m is imposed. If the cost C is expressed in terms of thermal energy units then m the function R( ) is simply a coefficient equal to 1.0. The cost of meeting the subsystem firm energy demand may alternatively be expressed in monetary units. This is particularly appropriate in the situation where the cost of thermal fuel (e.g., natural gas) varies considerably over the year. The cost function can be non-linear and time (i.e., month) dependent. However, when 33 used in conjunction with the formulation developed in Chapter 5, i t cannot vary from year to year. A multistage optimization problem whose objective is to minimize the cumulative (i.e., additive) costs arising at each stage and which has a single valued transformation function automatically meets the require-ments for decomposition and solution by Dynamic Programming [11]. If f (i ) is used to denote the cumulative future optimal cost mm c commencing at stage m and reservoir elevation im, then the set of recursive equations to be solved in the optimization of the hydrothermal subsystem operation over 12 monthly stages is as follows: f ^ i ^ = min[R(P1 - W(i ,d ,1 )j + (t ( i ^ d ,1^ ) ] I d l V V = mintR(Pm - " ' V W ' + fm + l ( t ( Vdm 'V ) ] 4 ' 5 . a m f 1 2 c i 1 2 > = min[R(F12 - w u 1 2 , d 1 2 f i 1 2 > ; + \iW\2Sr\l" The recursive solution of these equations is not possible until values are assigned to the optimal cumulative costs beyond the 12th stage, i.e., to the f ^ ( i ) for i = 1,. . . ,N. These values represent the terminal state values, discussed in Section 3.3, for operation terminating at the end of a current runoff year and reflect the influence of the reservoir state at the end of the year upon future operation. The costs assigned to the end state values are determined from the analysis of long-term operation formu-lated in Chapter 5. Equations (4.5) is solved by backwards recursion commencing with stage 12, determining the values of ^2 ^12^ f o r '''12 = * ', N' t n e n r ePe a t _ ing the process at stages 11, 10, 9, etc. until values of the optimal cumu-34 lative cost have been found for a l l reservoir states at stage 1. 4.2 RECOVERY OF OPTIMAL STATE TRAJECTORY AND ASSOCIATED COST Completion of the backwards (stage 12 to stage 1) recursive optimization phase of the Dynamic Programming analysis provides the opti-mal discharge decisions over a l l stages and for each state at the beginning of each stage. The sequence of monthly reservoir states which would occur i f the operation commenced at a particular state i ^ = i at the beginning of the year is obtained by applying the optimal decision for the i n i t i a l state and, using the state transformation Equation (4.2), determining the resulting state at the beginning of stage 2. This process is repeated over the remaining stages until the final state at the end of the last stage, i.e., the year end state, is obtained. If a l l aspects of the subsystem operation are regarded as fixed with the exception of the total annual inflow volume, which also defines monthly inflow volumes, then for a given i n i t i a l state i and an annual in-flow volume z, the associated year end state for the optimal trajectory is denoted as: j*(i,z) The cost incurred in operation of the subsystem on the optimal trajectory over the one year period is obtained by subtracting the optimal cumulative cost associated with the trajectory year end state from the optimal cumu-lative cost at the i n i t i a l state at the beginning of stage 1. Let the annual cost incurred along the optimal state trajectory for an i n i t i a l 3 5 state i and annual inflow volume z be denoted by: r*(i,z) then r*(i,z) = f ^ i ) - f 1 3 ( j * ( i , z ) ) ( 4 . 6 ) These optimal trajectory year end states and costs will be used in conjunc-tion with the analysis of long-term operation of the subsystem developed in Chapter 5. CHAPTER 5 FORMULATION OF THE ANNUL LONG-TERM OPERATION OF THE HYDRO THERMAL SUBSYSTEM The analytical scheme for the long-term operation of the hydro thermal subsystem is developed in this chapter by f i r s t showing that the sequence of annual operations, where the operating decisions depend upon the (random) annual inflow, constitutes a Markov decision process. It is then shown that the Dynamic Programming analysis of the within year operation described in Chapter 4 performs one phase of a two-phase iterative procedure which determines the optimal operating policy for the subsystem over an in-finite time horizon. 5.1 MARKOV PROCESS (OR CHAIN) IN DISCRETE TIME The basic feature of a discrete Markov process or chain (the Markov property) [17] is that the future behaviour of the process, which can be in one of a number of discrete states at successive discrete points in time, may depend upon the state of the system at a particular time but does not depend upon i t s earlier history. The discrete Markov process is specified by a matrix of one step transition probabilities and an i n i t i a l state (or a set of i n i t i a l state probabilities). The probability of a transition in one step from a state i at the be-ginning of a discrete time interval or stage to a state j at the end of the inter-val is a conditional probability: P[j|il 36 37 which is more conveniently written: P. . 13 The set of probabilities which define a l l possible outcomes from the one stage transition from an i n i t i a l state i may be written in the form of a row vector: { pi l Pi 2 Pi 3 * * ' Pi N} and the sum of the probabilities of a l l possible transitions from an i n i t i a l state must equal 1.0, i.e.: N .E. p.. = 1.0 3=1 13 The set of row vectors which defines the transition probabilities for a l l possible i n i t i a l to final state combinations during a single stage forms the one step transition probability matrix. Pl l * ' * Pl j ' * ' P1N o ft o • • • Pi l ' * ' Pi j * * ' Pi N « • « • • • PN1 ' * * PNj ' * * PNN A Markov process whose transition probability matrix does not change from one stage to the next is known as a Stationary or Homogeneous Markov process. These are usually the only Markov processes studied and the term Markov process as used throughout this thesis implies a stationary Markov process. 38 5.2 MARKOV PROCESS WITH RETURNS Suppose that a transition from one state to another involves a cost, a benefit, or more generally some form of quantifiable return being obtained. Let r.. denote the return for the one stage transition from an i n i t i a l state 13 i and a following state j . The complete set of transition returns for a l l possible i n i t i a l to final state transitions can be arranged in matrix form similar to the trans-ition probability matrix and is termed the transition return matrix: r i l ' * ' r i j ' * ' r i N • • • • • • r* r ' r* i l 13 lN « • • • • • • • • • • • 2-', _ • • • • r N 1 N3 NN 5.3 MARKOV DECISION PROCESS Decision making may be considered to be a process of selection of a particular action from a number of possible courses of action offering differing returns so that the return obtained is optimal in terms of some measurable criter-ion. Decision making is introduced into the Markov process by permitting a choice between a number of actions, each action resulting in different transition probabil-ity and return matrices. The process is further refined i f a decision consists of a set of actions with one action associated with each discrete state at the com-ic mencement of a transition. If d. is used to denote an action taken when state i 1 is encountered then the elements of the ith row of the transition probability and return matrices will be written in similar fashion: 39 k k Pi l Pi 2 k k ri l Pi 2 If there were K different action alternatives available at each of N the N states (i.e., k = 1,. . . ,K) then there would be K possible action com-binations and consequently KN associated transition probability and transition return matrices. As only stationary Markov processes are being considered here then only one of these action combinations is possible for the entire duration of any particular process. It is not necessary, however, that the same number or choice of actions be available for each state. k riN 5.4 DECISION PROCESS WITH RANDOM INPUT Consider a decision process with rewards where the state following is determined by a transformation function which is dependent upon the state i k prior to the transition, a decision d^, and an input Z which is an independent non-serially correlated random variable and takes on a discrete value z where z 1/. « . / B . k k Let t(i,d^,z) denote the transformation function and j (i,z) the state following this transition. Then: jk(i,z) = t(i,d*,z) 5.4.1 Similarly, let R(i,d^,z) denote the transition return function and r (i,z) the. actual return obtained from a transition from state i under decision 40 d. and random input z. Then: x rk(i,z) = R(i,d*,z) 5.4.2 The probability of a single transition from state i to state j follow-jc ing decision d^ is given by: p*. = P[j|i,djl = P[t(i,d*,z) = j|i,d*] 5.4.3 If the transformation function is inverted to the form z = t"1(i,d!C,j) 5.4.4 l and i f this inverted transformation uniquely defines a value of Z for given values of i , d^, and j (alternatively stated: i f only one random input can be jc responsible for a transition from i to j under decision d^ ) then the original and inverse equalities will both hold true coincidentally and therefore occur with the same probability. Equation 5.4.3 may therefore be rewritten: p* = P[z = t- 1( i , d * , j ) | i , d*] or Pi k = P [ Z = z | i' di] -1 k subject to; z = t ( i j d ^ j ) 5.4.5 but as Z is an independent random variable the probability of a particular value of Z occurring is not dependent upon the values of i and d^ and this equation further reduces to: 41 Pk. = P[Z=z] ID subject to: 5.4.6 z = t - ^ i . d ^ j ) Since the discrete probability distribution of the random input Z is known then the complete transition probability matrix is determined by Equation 5.4.6 for any specified decision set. Similarly, the transition return matrix is obtained from Equation 5.4.2. 5.5 DECISION PROCESS WITH INPUT DEPENDENT DECISIONS Now consider a decision process with random input as described above except that the decision enacted at the commencement of a state transition is dependent upon or related to the magnitude of the input occurring during that stage. The state transformation function must s t i l l uniquely define the following state as before so that once an input value has been revealed then only one state is feasible following the enactment of a particular decision. Because different decisions may now be made at an i n i t i a l state i (depending on the input value) then the same following state may now arise from two or more discrete input values and the inverse transformation function (Equation 5.4.4) no longer defines a unique value of z but now devines a set of values for z. k kz If D. represents the set of decisions {d. ; z = 1,. . . ,Z}, which l x contains the decisions associated with each of the discrete input values, then the set of discrete inflow values defined by the inverse transformation given the input dependent decision or policy set will be denoted by: {z = t'^i/D^j)} 42 The probability of a particular transition occurring is now the summation of the probabilities of a l l the possible events which will lead to that transition and Equation 5.4.6 now becomes; pk. = J ~ P[Z=z] 5.4.7 {z=t •L(i,D^,j)} Provided that the probability distribution of Z, the decision jc sets associated with a l l states D^ ; i = 1,. . .,N and the transformation function, do not vary for the duration of the process, then the one step transition matrix will remain the same for a l l transitions. Equation 5.4.7 indicates that the individual elements of the transition probability matrix are dependent upon the state of the system at a particular time but not upon its earlier history. The decision process with random input, where the de-cisions are dependent upon the value of the random input, therefore constitutes a stationary Markov process. A requirement of the transition probability matrix is that the sum of the probabilities of a l l possible transitions from any i n i t i a l state should be 1.0. That i s : _ k , n 5.4.8 E pi j = 1'° From Equation 5.4.7, - l k j {z=t (i,D ,j)} But: 43 L L j z=t-1<i,D*,j)} P [Z=z] p[z=z i : z j=t(i,Di,z) 5.4.10 and as the summation does not include a term which is a function of j then Ep* = E P[Z=z] 3 z = 1.0 5.5 LONG-TERM OPERATION OF THE HYDROTHERMAL SUBSYSTEM AS A MARKOV DECISION PROCESS The long-term operation of the hydrothermal subsystem may be con-sidered as a decision process with returns which undergoes a state transition during the course of each run-off year in the form of a change in reservoir storage volume. The annual inflow volume to the subsystem may be considered to be an independent non-serially correlated random variable comparable to the random input Z (see Section 3.2). This annual inflow is forecastable at the beginning of the runoff year so that a decision (in the form of an annual operating policy) which i s related to the forecast volume as well as the state of the reservoir at the beginning of the year can be prescribed, and enacted over the course of the year. Thus, provided a l l other aspects of the hydrothermal subsystem remain unchanged from year to year, the subsystem operation observed at yearly intervals has the elements of a Markov decision process with input dependent decisions and thus can be analyzed as such. In the remainder of this chapter, a procedure for the analysis of an unbound horizon Markov decision process is modified to suit the input dependent 44 process of the long-term operation of the hydrothermal subsystem. 5.6 DETERMINATION OF THE OPTIMAL LONG-TERM OPERATION POLICY FOR THE HYDROTHERMAL SUBSYSTEM Howard [18] has demonstrated a method for the determination of the optimal policy for an infinite horizon single recurrent chain Markov decision process using discounting. The method involves a two-phase iterative cycle consisting of a Policy Improvement routine and a Value Determination operation. In Section 5.6.1 the Dynamic Programming optimization of one year of subsystem operation, described in Chapter 4, will be shown to match the Policy Improvement routine when Howard's method is' applied to the long-term operation of the sub-system considered as a Markov decision process. The associated Value Determin-ation operation is then described in Section 5.6.2. 5.6.1 Policy Improvement Routine The Policy Improvement routine of the Howard method [18] evaluates the quantity; Hk = qk + S E pk. v(j) i r j 13 for each policy k and each discrete state i for i = 1,. . . ,N where: jc q. = expected return from a l l possible one stage transitions from state i ; 8 = discount factor for the duration of one stage; v(j) = the present worth of the total expected return for the process in state j , j = 1,. . .,N * The policy k is then selected as the one which, in the case of a mini-k * mization problem, produces the lowest value of H., namely H.. Thus: 45 * k k H. = minimum [q. + g£p..v(j)] 5.6.1 i i j ID For the case of the Markov decision process with discrete random input Z which takes on value z;z = 1,. „ .,3, and decision dependent upon this random input as discussed in Section 5.5., i t was shown that: Given a state transformation jk(i,z) = t(i,Dk,z) and one stage transition return rk(i,z) = R(i,Dk,z) then the transition probability for this process was given by Equation 5.4.7 1 3 { z - t : < i , D ? , j ) > 5.6.2 and the expected one stage return commencing at state i will be the summation of the products of the transition probabilities and their associated returns: =1 Y. _! k P[Z=z]rk(i,z) j {z=t (ijD^j)} 5.6.3 k k Substituting for p ^ and q^ from Equations 5.6.2 and 5.6.3, then Equation 5.6.1 may now be written: H. = minimum 1 *k i j {z=t (i,D.,j)} P[Z=z]rk(i,z) P[Z=z]v(j) j {z=t (i,Di,j)} 46 = minimum 1 I I P[Z=z]{r (i,z) + 6v(j)} j { z ^ - ^ i , ^ , j ) } 5.6.4 The double summation may be replaced by an equivalent double sum-mation (see Equation 5.4.10) so that Equation 5.6.4 becomes: H. = minimum l P[Z=z]{rk(i,z) + 8 v(j)} z j=t(i,Difz) As the transformation equation in its non-inverted form is single k k valued and j = t(i,D^/z) = j (i,z) then: H. = minimum 1 ^>P[Z=z] (rk(i,z) + 8 v(jk(i,z))} 5.6.5 As the quantity within the parenthesis {} is independent of the probabil-k kz ity of Z=z and as represents a set of independent decision sets each associated with a state and inflow combination then the minimization operation may be brought within the summation and confined to selecting the optimal de-cision for a particular input and state combination. Hence: * V k v H. =2_ptz=zl minimum {r (i,z) + 8 v(j (i,z))} X z j c z 5.6.6 Thus each minimization is performed for a particular state i , dis-crete input z, and set of state values v(j) for j =1,. . . ,N and minimizes the sum of the cost over a single transition and the discounted value of the state following the transition. In the context of the long-term operation of the hydrothermal subsystem this is equivalent to determining the minimum cumulative cost of operation (i.e., 47 a n n u a l o p e r a t i n g c o s t p l u s t h e p r e s e n t w o r t h o f t h e e x p e c t e d c o s t o f f u t u r e o p e r a t i o n ) f o r a g i v e n r e s e r v o i r s t a t e i a t t h e b e g i n n i n g o f a y e a r a n d a n a n n u a l i n f l o w v o l u m e f o r t h e y e a r z . T h i s m i n i m i z a t i o n o b j e c t i v e i s i d e n t i -c a l t o t h a t u s e d i n t h e o p t i m i z a t i o n o f o n e y e a r o p e r a t i o n o f t h e s u b s y s t e m i n C h a p t e r 4 w i t h t h e e x c e p t i o n t h a t i n t h e l a t t e r t h e a n n u a l c o s t w a s d e -c o m p o s e d i n t o t h e s u m o f t h e 12 m o n t h l y c o s t s a n d t h e y e a r e n d s t a t e v a l u e s w e r e n o t d i s c o u n t e d . U n d e r t h e c o n d i t i o n t h a t t h e m o n t h l y i n f l o w s a r e u n i q u e l y d e f i n e d b y t h e a n n u a l i n f l o w v o l u m e a n d r e p l a c i n g t h e t e r m i n a l s t a t e v a l u e s s o t h a t : F 1 3 ( j ) = $ v ( j ) f o r j = 1,. . . , N t h e n t h e D y n a m i c P r o g r a m m i n g a n a l y s i s o f t h e o n e y e a r o p e r a t i o n a s f o r m u l a t e d i n C h a p t e r 4 a n d p e r f o r m e d f o r a l l d i s c r e t e i n f l o w v a l u e s w i l l p e r f o r m t h e o p t i m i z a t i o n r e q u i r e d f o r t h e P o l i c y I m p r o v e m e n t r o u t i n e i n t h e a n a l y s i s o f l o n g - t e r m s u b s y s t e m o p e r a t i o n . I f r * ( i , z ) a n d j * ( i , z ) a r e a g a i n u s e d t o d e n o t e t h e a n n u a l c o s t a n d y e a r e n d s t a t e a s s o c i a t e d w i t h t h e o p t i m a l s t a t e t r a j e c t o r y f o r i a n d z o b t a i n e d f r o m t h e m o d i f i e d D y n a m i c P r o g r a m m i n g a n a l y s i s t h e n E q u a t i o n 5.6.6 m a y b e w r i t t e n : H ? = ^ P [ Z = z ] { r * ( i , z ) + B v ( j * ( i , z ) ) } 5.6.7 z T h e r e s u l t s o f t h e o p t i m i z a t i o n a l s o i n c l u d e t h e o p t i m a l a n n u a l d e c i s i o n s e t s f o r e a c h s t a t e , l e t t h e s e b e d e n o t e d b y D ^ , f o r i = 1,. . . , N . * T h e e l e m e n t s p . . o f t h e m a t r i x o f o n e y e a r t r a n s i t i o n p r o b a b i l i t i e s a s s o c i a t e d 48 with are obtained from the probability distribution of Z and Equation 5.6.2: Ji i - = 2- i * 1 3 {z=t-1(i,D.,j)} P[Z=z] for i = 1,. . .,N; j = 1,. . .N Similarly, the expected one year transition costs q^ associated * with are obtained from the probability distribution of Z, the values of * r (i,z) and Equation 5.6.3: " E X . ! * P[Z=z]r* j (z=t X(i,D.,j)} (i,z) for i = 1,. , .,N 5.6.8 The resulting one year transition probability matrix and vector of expected one year transition costs are utilized in the Value Determination operation described in the following section. 5.6.2 Value Determination Operation The Value Determination phase of the Howard method computes the present worth expected cost v'(j) for each state j = 1,. . ,,N for the Markov process having the transition probability matrix and vector of expected one stage transition returns obtained from the Policy Improvement routine described in Section 5.6.1. These state values are obtained by solving the set of simultaneous linear equations: v'(i) = q. +8 E p.. V ( j ) for i = 1,. . „,N 5.6.9 3 The resulting state values are then returned to the Policy Improvement routine where the Dynamic Programming analysis of one year operation of the hydrothermal subsystem is repeated with these state values. A revised trans-4 9 ition probability matrix and vector of one stage transition returns are calculated and the Value Determination operation described above is re-peated. This iterative process is continued until the state values ob-tained by the Value Determination operation on two successive iterations are the same. When this occurs the operating policy cannot be further im-proved and is optimal for the long-term operation. A proof of the conver-gence properties of the iterative cycle is given by Howard [18] and is un-affected by the modified Policy Improvement routine described in this chapter. 5.6.3 Long-Term Operating Cost With Unspecified Initial State The state value v(j) provides a value for the present worth expected cost of long-term operation of the hydrothermal subsystem commencing at a specified i n i t i a l reservoir state j . If a value of the cost of long-term operation is required when an i n i t i a l reservoir state cannot be specified (e.g., in the case of an incomplete project) then a representative cost may be determined from the sum of the state values for the recurrent states weighted by their associated steady state probabilities, i.e., the sum: E 1L • v(j) 5.6.10 where 11.. is the steady state probability for state j determined from the opti-mal transition probability matrix | p „ | . The result of this summation will be referred to as the steady state present worth expected cost. CHAPTER 6 COMPUTER PROGRAM AND COMPUTATIONAL EXPERIENCE The computational scheme for determining the optimal long-term operation of the hydrothermal subsystem for a specific annual firm energy output and thermal capacity allocation follows the two phase iterative cycle described in Sections 5 . 5 and 5 . 6 . The program consistsof a main program which includes the Value Determination stage of the cycle, a Dynamic Programming sub-routine which performs the optimization of the Policy Improvement phase, and a simulation subroutine which determines the results of operating the hydroelec-tric project over a one-month period given a particular decision and the oper-ating conditions. A general description of the main computer program and subroutines, with computer flowcharts, together with a summary of computational experience with the program, is given in this chapter. A complete listing of the computer program is given in Appendix 1. 6 . 1 MAIN PROGRAM A flow chart for the main program is given in Figure 6 . 1 . Initial values of the following key variables are specified for each computation of optimal long-term operation: Annual discount factor. Upper and lower limits of reservoir level. Required firm energy output from the hydrothermal subsystem. An i n i t i a l set of state values v(j) for j = 1 , . . .,N. 50 (START ) Read Subsystem Data 51 T POLICY I M P R O V E M E N T Dynamic Programming Optimalisation ( Subrout ine OP ) I V A L U E D E T E R M I N A T I O N A. Solve L i n e a r Eqns . for v ( j ) Has Conversance" been j ichievec Y e s No Ca lcu la te S teady State Probab i l i t i es Ca l cu la te Present Worth Expected Cost Print results for this Firm Energy C J u t p u L ^ ^ • Increment Firm Energy Output (STOP ) FIGURE 6 .1 52 The following basic data are also provided; Values of discrete annual inflows and their associated probabilities of occurrence. Reservoir elevation and storage relationships. Coefficients for determining monthly inflows from annual inflow. Coefficients for shaping the monthly energy demand. The iterative cycle commences by calling the DP subroutine which returns a vector of expected transition returns and the transition probability matrix given the i n i t i a l set of state values, v(j); j = 1,. . .,N. This infor-mation is then used in the value determination portion of the main program to determine a new set of state values. A comparison is made between the previous set of state values and the newly computed ones; i f these do not agree to within a specified tolerance (i.e., i f convergence of the state values has not been achieved at that iteration), the program returns the latest set of state values to the DP subroutine and the cycle continues. When convergence of the state values occurs then the steady state probability vector is obtained by progressively raising the power of the final transition probability matrix obtained from the DP subroutine until the result-ing matrix stabilizes, i.e., when: The long-term present worth expected cost for steady state operation at the firm energy level and thermal plant capacity prescribed in the i n i t i a l values is finally computed by multiplying the optimal steady state probability vector by the vector of optimal state values. 53 6.2 DP SUBROUTINE The DP subroutine determines the optimal twelve month operating trajectories for a l l i n i t i a l states and for a l l discrete annual inflows. The computer flowchart is given in Figure 6.2. The optimization procedure for a particular annual inflow is performed in two phases: determination of the optimal decision for a l l states and stages from the backward recursive equations 4.5 and the current set of state values v(j); j = 1,. . ,,N, and then retracing the optimal trajectory for a l l i n i t i a l reservoir state values. The optimization involves the determination of the outcome of each feasible decision at each state and stage combination in terms of hydro power generation, thermal energy requirement (hence cost) and resulting following reservoir state. This is performed in the Simulation subroutine. The optimal trajectories are recovered by tracing the sequence of optimal decisions and following states commencing with a state at the begin-ning of the firs t stage. Operating costs are accumulated at each stage along the trajectory and the total over the 12 stages constitutes the one year opti-mal transition cost. The reservoir state at the termination of the one year optimal trajectory is the optimal terminal state for the transition from the state at the beginning of the year. The transition probability matrix and the vector of expected transition costs are determined from this trajectory infor-mation and the annual inflow probabilities. (START) In f low Loop Stage Loop S ta te Loop Dec i s i on Loop I Simu la te e f f e c t of D e c i s i o n at this s ta te a s t a g e ( S B R T N . S I M ) No Update Optimal Decision No No No 54 ® FIGURE 6.2 COMPUTER FLOWCHART FOR DYNAMIC PROGRAMMING SUBROUTINE (Continued on Following Page) O P T I M A L T R A J E C T O R Y R E C O V E R Y Stage Loop In i t ia l S ta te Loop I R e c o v e r Opt ima l cos t , State f o l l o w i n g , Opt imal d e c i s i o n for Opt imal Tra jec tory c o m m e n c i n g at I n i t i a l S t a t e No No No Dete rm ine E x p e c t e d one Y e a r T r a n s i t i o n c o s t s , P r o b a b i l i t i e s X (RETURN) FIGURE 6.2 (Continued) 56 6.3 SIMULATION SUBROUTINE See computer flowchart Figure 6.3. The simulation subroutine is provided with the following infor-mation from the main and DP (calling) programs-Stage (month) and reservoir state. A decision value, i.e., the level(state) at which the reservoir is to be at the end of the current stage. Firm energy demand for the month. Inflow for the month. State-reservoir storage relationships. The decision is defined within the DP subroutine as a state change over the month stage. The simulation subroutine computes the discharge which would be necessary to produce the state change defined by the decision, allow-ing for the reservoir inflow volume during the month. The program then com-putes the amount of thermal energy which would be required to augment the hydro output and meet the energy demanded. A check is made to ensure that a l l dis-charges, energy outputs, and reservoir elevations are feasible with respect to the relevant capacity and reservoir constraints. If any constraint is violated then an extremely high cost (e.g., 10*° thermal units) for that decision is re-turned to the DP subroutine. The magnitude of the penalty is made large enough to ensure that an optimal trajectory which includes an infeasible decision can be recognized in the main program. 57 ( START) Dete rm ine D i s c h a r g e V o l u m e R e q u i r e d C a l c u l a t e Hydro Output No C a l c u l a t e The rma l E n e r g y R e q u i r e d to meet Firm Energy Demand No C a l c u l a t e cost of thermal ene rgy ( o p t i o n a l ) (RETURN) A l l o c a t e Penalty C o s t = I O l 2 u n i t s -< F I G U R E 6 . 3 C O M P U T E R F L O W C H A R T F O R S I M U L A T I O N S U B R O U T I N E 58 6.4 COMPUTATIONAL EXPERIENCE 6.4.1 Convergence of the Iterative Cycle In order to achieve greater sensitivity in the selection of the optimal monthly discharge decisions the computer program was originally written so that decisions which resulted in a transition to an intermediate non-discrete state could be considered. Solution of the recursive optimization involved the use of an interpolative procedure to determine the value of the optimal future cost for the intermediate state at the beginning of the next stage. Even with simple linear interpolation procedures this has been shown to provide much greater accuracy than discrete Dynamic Programming [24] . However, when i t is necessary to recover optimal trajectories, as was the case here, then further interpolation for optimal decision values at intermediate states is necessary. Unfortunately i t was found that with this double interpolation the iterative cycle had very unstable characteristics and showed no tendency to converge on an optimal solution even after a large number of iterations. A series of damp-ing devices were introduced in an attempt to suppress the instability but even so convergence could not be consistently achieved. Consideration of the effect of small errors in the iterative cycle suggests why this problem arose. If we consider the point in the iterative cycle where the policy improvement phase has returned a truly optimal vector of ** ** expected transition returns Q and transition probability matrix P . The solution of the linear equations in the value determination phase will then ** give an optimal set of endstate values V but will introduce small errors, ** Ihus i f V represents the vector of true state values for the process described ** ** ** by Q , P and e the vector of errors, then the actual vector of state values computed will be; 59 ** ** V + e In order to confirm that convergence has been achieved the iterative cycle must be repeated once more. Errors in the state values will in turn produce errors in the expected transition returns and transition probability matrix computed in the subsequent policy improvement phase, let these be represented by: ** ** Q + e' and P + e" The set of linear equations which are then solved to produce the new set of state values for this cycle will now be: V'(i) = (q** + £i) + 6E(p** + e^ .) V (j) i = 1,. . .,N When solved, the v' ( i ) , which represent the present worth of the expected cost of operation over an unlimited period of time, will be large compared with the q^ term, which represents the expected cost of operation over a single year, i and small additive error in the expected transition returns will not induce greater order errors. In the summation term, however, i t can be seen that errors occurring in elements of the transition probability matrix are multiplied by the large magnitude state values and will result in magnified errors in the v'( i ) . Since the use of interpolative DP results in year-end states which do not necessarily coincide with discrete state values these must be made dis-crete in some way. A slight shift in an optimal trajectory due perhaps to a small error in the iterative cycle may result in the assignment of a year-end state from one discrete state value to another, i.e., induce a change equal to 60 the interval between discrete states. While this may have negligible effect on the expected transition cost for the annual trajectory i t will cause a significant change in the probability transition matrix. In view of the sensitivity to error in the transition probability matrix as demonstrated above state interpolation was considered to be the prime cause for failure to converge. The Dynamic Programming subroutine was rewritten in fully discrete form with no interpolation. The decision variable being in the form of a decision to change the state of the reservoir an integer number of state in-crements thus ensuring that state transformations would coincide with discrete following states. This revision immediately provided the reproducibility re-quired to achieve consistent convergence but at the cost of some loss of accur-acy in the analysis. 6.4.2 Initialization and Computational Time The set of state values v(j) for j = 1,. . .,N used to initialize the iterative cycle may be chosen quite arbitrarily. The DP subroutine can only select feasible state trajectories and hence must produce reasonable trans-ition costs. Consequently the fi r s t value determination computation will always provide a realistic set of state values. It was found that approximately six iterations were required to achieve convergence when a purely arbitrary set of i n i t i a l state values was used. When the computation was arranged so that the optimal state values obtained for a lower or higher firm energy level were used as i n i t i a l state values, then i t was found that the number of iterations required in the current computation was halved or usually about three. In order 61 to obtain one of the principal results in this study (see Figure 8.1) optimi-zation of the hydrothermal subsystem is required at progressively increasing levels of annual firm energy demand. It was therefore only necessary to pro-vide arbitrary initializing state values for the f i r s t firm energy level compu-tation and i t was then possible to initialize subsequent computations with the values obtained from the previous (lower) firm energy optimal solution. This 50 per cent reduction in computational time achieved when producing the results shown in Figure 8.1.1 is significant since i t permits the use of smaller state and decision increments and hence greater accuracy for acceptable compu-tational times. Computational time on the IBM 360/67 for the 20 discrete state, 12 stage and 9 discrete annual inflow analysis was approximately 10 seconds per iteration or 30 seconds per point when plotting the firm energy/thermal cost graph referred to above. CHAPTER 7 STATISTICAL DATA FOR THE PORTAGE MOUNTAIN SUBSYSTEM 7.1 TEST FOR SERIAL DEPENDENCE OF ANNUAL INFLOWS A requirement of the analysis presented in Chapter 5 is that the annual inflow to the overyear storage reservoir mus tbe an independent non-serially correlated random variable. The historic record of 39 annual in-flows was tested for serial dependence by fi r s t fitting (least squares re-gression) a f i r s t order serial correlation model relating the inflow in year I with the inflow in the previous year I ,. This model is of the form: y y-i I = bn + b.I + e 7.1.: y 0 1 y-1 The results obtained for 38 years pairs 1,1 , obtained from the 39 years y y-1 of historic record yielded a f i r s t order serial correlation coefficient r^ of 0.2361. The test of significance developed by Anderson [19] based on a ran-dom normal time series gives confidence limits for the value of the serial correlation coefficient. At the 95 per cent confidence level for this sample size the confidence limits for the f i r s t order serial correlation coefficient r^ are -0.345 to +0.291. Since the computed value of r^ falls between these limits i t is considered to be insignificantly different from zero at the 95 per cent confidence level and therefore suggests that no significant corre-lation exists. 62 63 A second test was applied in the form of the Durbin-Watson d statistic test for serial correlation of the error terms ofthe least squares regression model Equation 7.1.1 [23]. The computed value of the Durbin-Watson d statistic was 2.126 which lies above the upper limit of du = 1.54 for this sample size at the 95 per cent confidence level. This result also confirms that there is no significant serial correlation in the annual in-flow volumes. As there was no other evidence to suggest that significant serial correlation exists between annual inflow volumes to the Portage Mountain project i t was concluded that the annual inflow could be reasonably assumed to be an independent non-serially correlated random variable. 7.2 PROBABILITY DISTRIBUTION OF THE ANNUAL INFLOW VOLUME In the absence of any prior statistical knowledge which might pro-vide a theoretical justification for a particular probability distribution for the annual inflow volume a curve fitting procedure was chosen. The Pear-son system of frequency curves was selected to provide the mathematical model of the distribution. Although this model has only slight theoretical basis [19] i t has been widely used and provides a systematic and straightforward procedure for fitting a four parameter curve to virtually any distribution. The fitting procedure [20] commences with the calculation of the fi r s t four moments of the sample data and then computes the value of a c r i -terion k which determines the type of Pearson curve to be fitted. The value 64 of k obtained from the annual inflow data indicated a Pearson Type I d i s t r i -bution. This is a bell shaped curve, skewed, and with limited range in both directions. The values obtained for the various parameters and constants in-volved in the fitting procedure are given in Table 7.2.1 and the resulting frequency distribution curve is shown in Figure 7.2.1. 2 The fitted curve was tested for goodness of f i t by both the x and the Kolmogorov-Smirnov test. The latter test is considered the more power-ful when the hypothesized distribution is fully defined and, by comparing cumulative distribution, avoids the necessity of grouping data into class intervals. The results obtained from both tests confirmed that the fitted curve agreed with the observed data at the 99 per cent confidence level. Nine inflow classes were selected over the finite range of the fitted distribution curve Figure 7.2.1. The probability of an inflow occurring in each of the inflow classes was then determined together with the weighted average inflow for each class. The inflow class intervals class probabilities, and weighted average inflow for each class are given in Table 7.2.2. These results are used to provide the discrete input values and their associated probabilities in the suboptimization analysis described in Chapter 5. FIGURE 7.2.1 FREQUENCY DISTRIBUTION FOR ANNUAL INFLOW VOLUME PORTAGE MOUNTAIN SUBSYSTEM 66 TABLE 7.2.1 RESULTS OF FITTING PEARSON TYPE I CURVE TO PORTAGE MOUNTAIN HISTORIC ANNUAL INFLOW VOLUMES FIRST FOUR MOMENTS FROM DATA: 1st Moment 2nd Moment 3rd Moment 4th Moment 1.147 0.0179 0.892 x 10 0.029 x 10 RESULTS OF FITTING PROCEDURE: Criterion -0.173 (negative confirms Pearson Type I curve) Mode 1.116 Intercepts Left 0.769 Right 1.998 EQUATION OF FITTED PEARSON TYPE I DISTRIBUTION CURVE: /« ,i.x.m_ x %m p(x) -p o(l +-) 1 d - - ) 2 f o r _ a i < x < a 2 p ( x, = 2 . 9 1 ( i + _ | - _ ) ^ 2 3 ( 1 . _ x _ _ ) 2 . 9 1 where x is measured from origin at mode. 12 Note: Where applicable units are 10 cu.ft. 67 TABLE 7.2.2 DISCRETE ANNUAL INFLOW VOLUMES AND PROBABILITIES ANNUAL INFLOW VOLUME CLASS LIMITS Left 10 cu.ft. Right CLASS WEIGHTED 12" Average 10 cu. f t . Probability 0.769 0.900 0.876 .016 0.900 1.000 0.960 .121 1.000 1.100 1.053 .254 1.100 1.200 1.148 .278 1.200 1.300 1.245 .196 1.300 1.400 1.432 .095 1.400 1.500 1.439 .032 1.500 1.600 1.535 .007 1.600 1.998 1.635 .001 1.000 68 7.3 DETERMINATION OF MONTHLY INFLOWS FROM A TOTAL ANNUAL INFLOW In order to develop a deterministic relationship between the total annual inflow volume and the monthly inflow values which make up this annual volume a linear regression model was adopted in the form of Equation 7.3.1: X = b__ + b, (Y) for m = 1,. . .,12 7.3.1 m Om lm where m denotes a particular month, b„ and b, are the regression coeffi-Om lm cients for estimating the inflow volume X in month m from the annual in-m flow volume Y. The 39 year record of monthly inflows was used to provide 39 sample pairs of X^ and Y for fitting the regression equation for each month. The values obtained for the regression coefficients and correlation coeffi-cients for each month are given in Table 7.3.1. Significant correlation was obtained at the one per cent level for the months of June, July, August and September. As the average combined in-flow over these months constitutes on the average approximately 70 per cent of the total annual inflow, this monthly inflow model was considered as sat-isfactory. A requirement of the set of 12 monthly inflow equations is that for a given annual inflow volume Y they will generate a set of monthly in-flows whose total is equal to Y. The proof that this requirement is met is given in Appendix 2. 69 TABLE 7.3.1 REGRESSION AND CORRELATION COEFFICIENTS FOR GENERATION OF MONTHLY INFLOW VOLUMES REGRESSION COEFFICIENTS Om lm CORRELATION COEFFICIENT r January February March April May June July August September October November December 0.3 1.2 5.6 38.5 164.2 -80.0 -32.2 -60.7 -17.1 - 0.4 -19.1 - 0.4 0.0186 0.0170 0.0119 0.0058 0.0379 0.3419 0.1903 0.1423 0.0781 0.0647 0.0637 0.0280 0.37 0.28 0.27 0.04 0.10 0.66 0.62 0.65 0.70 0.48 0.50 0.37 70 7.4 MONTHLY LOAD DISTRIBUTION The distribution of the annual firm energy demand by month which was imposed upon the Portage Mountain subsystem was based upon the firm energy forecast data for the overall, generating system [21]. The distribution was adjusted for uniform month length (i.e., 730 hour months) and is given in Table 7.4.1. Also shown in Table 7.4.1 is the monthly distribution of the com-bined output (i.e., proportion of annual output generated in each month) of the existing seasonal storage and run of river hydro projects in the B.C. Hydro generating system. This data has been obtained from computer simula-tion of the overall generating system [21]. This data was not analyzed but from inspection exhibits apparent random fluctuation both from month to month and from year to year which suggests that the effect on the residual firm energy load imposed on the overyear storage hydro project(s) is to intro-duce a noise component but not a consistent pattern which is useful in predic-ting load shaping. The annual firm energy load shape suggested by the fore-cast data was therefore adopted in this study. TABLE 7.4.1 DISTRIBUTION OF PROPORTION OF ANNUAL FIRM ENERGY OUTPUT BY MONTH MONTH PROPORTION OF ANNUAL FIRM ENERGY OUTPUT FOR ALL HYDRO PROJECTS EXCLUDING OVERYEAR STORAGE HYDRO* PROPORTION OF ANNUAL SYSTEM FIRM ENERGY DEMAND FROM FORECAST DATA January February March April May June July August September October November December Average .09 .07 .08 .07 .09 .09 .09 .08 .07 .08 .09 .09 Max. .11 .09 .09 .09 .10 .11 .11 .09 .09 .10 .11 .11 Min. .07 .06 .07 .06 .07 .08 .07 .07 .06 .06 .07 .07 .091 .094 .090 .079 .076 .076 .074 .077 .079 .082 .091 .091 * Obtained from a computer simulation of system [21]„ CHAPTER 8 RESULTS OF THE ANALYSIS OF THE LONG-TERM OPERATION OF THE OVERYEAR STORAGE SUBSYSTEM The results presented in this chapter were obtained from the analysis of the long-term operation of the hydrothermal subsystem which included the Portage Mountain project of the BCHPA generating system. 8.1 ANALYSIS OF THE HYDROTHERMAL SUBSYSTEM WITHOUT THERMAL CAPACITY CONSTRAINT A single run of the computer program described in Chapter 6 produced, for a prescribed annual firm energy output F, an optimal set of state values v(j) for j = 1,. . .,N and the present worth expected cost (PWEC) of steady state operation. If the elevation of the storage reservoir at the beginning of the current runoff year is known then the state value for the state corresponding to that elevation is the most appropriate measure of the PWEC for the subsystem operation at that firm energy output. Alternatively, i f for example, the project is not yet operational, then the PWEC for steady state operation computed according to Equation 5.6.10 may be used. It should be noted that, for the purposes of this study, the PWEC is measured in terms of thermal energy units consumed and is not converted to monetary units as discussed in Section 4.1. The results obtained when the analysis of long-term operation (i.e., the suboptimization process) is performed over a range of firm energy outputs for the subsystem and with no constraint on the thermal energy required in 72 7 3 any month are summarized in Tables 8.1.1 to 8.1.6. These results are more conveniently shown in graphical form in Figure 8.1.1 which was ob-tained by plotting the subsystem annual firm energy output against the present worth expected cost for steady state operation. Similar graphs are obtained i f the state value for a particular reservoir state is plotted at each firm energy output and examples of these are also shown in Figure 8.1.1. For the purposes of the discussion following the PWEC for steady state operation will be used. The annual firm energy output F.^ on Figure 8.1.1 represents the maximum firm energy output which can be achieved by the hydrothermal sub-system without consumption of any thermal energy and therefore is equivalent to non-integrated operation of the hydro project alone. The line ac repre-sents the cost of integrated operation for subsystem firm energy outputs above F^ and is the principal result obtained from the suboptimization. The line ab represents the subsystem performance i f the operation of the hydro and thermal plants within the subsystem is not integrated and is obtained from a simple energy balance equation. If F represents an annual firm energy output and T represents the annual consumption of thermal energy by the subsystem when the operation of the thermal and hydro plant is not integrated, then for F >_ F^: F = F + T 8.1.1 or T = F - F. 8.1.2 74 TABLE 8.1.1 SUBOPTIMIZATION RESULTS• PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 10000. mkwh STEADY STATE RESERVOIR STATE STEADY STATE PRESENT WORTH STATE VALUE PROBABILITY EXPECTED COST j V ( j ) 7 T ( j ) C 20 (el. 2200 ft.) 30.7 19 30.8 .008 18 31.0 .116 17 31.5 .303 16 33.0 .055 15 35.6 .158 14 39.3 .080 13 43.7 .082 12 50.7 .056 11 59.2 .043 10 69.1 .029 9 78.9 .031 8 108.1 .038 7 396.3 6 720.1 5 1079.3 4 1441.3 3 1802.2 2 2156.1 1 (el. 2150 ft.) 2502.9 1.000 TABLE 8.1.2 SUBOPTIMIZATION RESULTS: PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 12000. mkwh STEADY STATE RESERVOIR STATE STEADY STATE PRESENT V70RTH STATE VALUE PROBABILITY EXPECTED COST j v(j) ir(j) C 20 (el. 2200 ft.) 3935.5 19 3939.6 18 3964.2 0.266 17 4024.2 0.208 16 4147.6 0.123 15 4286.3 0.143 14 4457.2 . 0.149 4365.7 13 4639.5 0.019 12 4842.3 0.025 11 5044.3 0.040 10 5250.1 0.015 9 5473.2 0.012 8 5708.6 7 6092.7 6 6476.6 5 6853.6 4 7222.7 3 7585.0 2 7939.5 1 (el. 2150 ft.) 8287.9 1.000 TABLE 8.1.3 SUBOPTIMIZATION RESULTS: PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 14000. mkwh STEADY STATE RESERVOIR STATE STEADY STATE PRESENT WORTH STATE VALUE PROBABILITY EXPECTED COST j v ( j ) ir(j) C 20 (el 2200 ft.) 18000. 19 18117. 18 18296. 17 18510. .226 16 18789. .171 15 19114. .094 14 19446. .133 13 19774. .128 12 20100. .024 11 20424. .086 10 20763. .065 9 21130. .036 8 21506. .036 7 21890. 6 22274. 5 22651. 4 23021. 3 23383. 2 23737. 1 (el. 2150 ft.) 24086. 1.000 TABLE 8.1.4 SUBOPTIMIZATION RESULTS: PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 16000. mkwh RESERVOIR STATE STATE VALUE STEADY STATE PROBABILITY STEADY STATE PRESENT WORTH EXPECTED COST j v(j) TT ( j ) 20 (el. 2200 ft.) 38572.9 19 38849.1 18 39163.0 17 39517.8 0.217 16 39933.2 0.204 15 40390.7 0.442 14 40830.2 0.121 13 41273.6 0.016 12 41710.9 11 42140.1 10 42561.4 9 42974.7 8 43379.2 7 43768.0 6 44156.7 5 44539.6 4 44916.6 3 45289.2 2 45654.8 1 (el. 2150 ft.) 46012.3 40174.8 1.000 78 TABLE 8.1.5 SUBOPTIMIZATION RESULTS: PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 18000. mkwh STEADY STATE RESERVOIR STATE STEADY STATE PRESENT WORTH STATE VALUE PROBABILITY EXPECTED COST j V ( j ) TT ( j ) C 20 (el. 2200 ft.) 63311.3 19 63751.2 18 64205.6 17 64666.4 0.040 16 65125.4 0.095 15 65612.4 0.135 66044.4 14 66068.4 0.265 13 66517.0 0.228 12 66957.2 0.177 11 67388.6 10 67811.8 9 68227.7 8 68637.0 7 69030.2 6 69428.2 5 69821.3 4 70207.4 3 70586.3 2 70957.6 1 (el. 2150 ft.) 71323.0 1.000 79 TABLE 8,1.6 SUBOPTIMIZATION RESULTS: PORTAGE MOUNTAIN SUBSYSTEM SUBSYSTEM ANNUAL FIRM ENERGY OUTPUT 20000. mkwh STEADY STATE RESERVOIR STATE STEADY STATE PRESENT WORTH STATE VALUE PROBABILITY EXPECTED COST j v(j) ir(j) C 20 (el. 2200 ft . 89291.1 19 89773.4 18 90253.4 17 90731.8 0.001 16 91205.1 0.007 15 91696.0 0.080 14 92154.6 0.277 13 92603.9 0.312 12 93045.4 0.168 11 93479.5 0.102 10 93904.6 0.035 9 94322.4 0.018 8 94735.3 7 95132.9 6 95535.2 5 95933.1 4 96323.8 3 96707.6 2 97082.6 1 (el. 2150 ft.) 97450.8 1.000 Present Worth Expected Cost mkwh (thermal energy) o O 81 As represents the firm energy capability of the hydro project under any (discrete) annual inflow situation its value is constant for a l l years. When the annual consumption of thermal energy by the subsystem T is plotted against the annual subsystem firm energy output F then from Equation 8.1.1 the slope of the graph for non-integrated operation will be a line of slope 1:1. The present worth expected cost for non-integrated operation will therefore be the present worth of a uniform annual series for an amount T per year (with probability 1.0) commencing at the beginning of the f i r s t year and extending indefinitely into the future. The present worth factor for such a series corresponding to the annual discount factor 0.926 used in the analysis is 13.5. If the present worth expected cost of non-intregrated. operation with an annual firm energy output F is denoted by PWEC then: F PWEC = 13.5 T F Hence from Equation 8,1.2 PWEC = 13.5 (F - F.) F 1 and the slope of the line ab in Figure 8.1.1 will be 13.5:1. The vertical distance between the line ab and the curve ac at an annual firm energy output is a measure of the reduction in cost, specifi-cally the present worth expected cost, achieved from integrated operation. This cost reduction increases progressively for increasing values of F beyond F^ until point d is reached beyond which the cost reduction remains constant and the integrated and non-integrated lines become parallel. This indicates 82 that beyond the annual firm energy output corresponding to point d any additional firm energy output demanded from the subsystem can only be ob-tained at the same incremental cost as would be incurred in the non-integrated operation of a thermal plant. 8.2 ANALYSIS OF THE HYDROTHERMAL SUBSYSTEM WITH THERMAL GENERATING CONSTRAINT INCLUDED When a thermal capacity constraint in the form of a maximum thermal energy generation in the subsystem in any month is imposed then its perform-ance is modified. The results obtained for progressively increasing values of F with three different values of the thermal capacity constraint are given in Table 8.2.1 and are shown graphically in Figure 8.2.1. The uncon-strained or infinite thermal capacity curve and the non-integrated operation line from Figure 8.1.1 are also superimposed on this graph. The integrated performance curve for a particular thermal capacity T mkwh/month (lmkwh = 106 kwh) follows the infinite capacity curve until the annual firm energy out-put approaches the maximum which can be achieved by the subsystem F^ regard-less of cost (for that thermal capacity T). In this region the cost rises rapidly and the curve terminates when i t intersects the non-integrated line ab at an annual firm energy output F . If T^ represents the maximum thermal energy which can be generated in a year with a thermal capacity of x mkwh/month so that: T = 12T 8.2.1 T then the value of F is obtained from Equation 8.1.1: TABLE 8.2.1 SUBOPTIMIZATION RESULTS: PORTAGE MOUNTAIN SUBSYSTEM WITH FINITE THERMAL CAPACITY THERMAL CAPACITY mkwh/month ANNUAL FIRM ENERGY OUTPUT mkwh STEADY STATE PRESENT WORTH EXPECTED COST 9500 4.3 9750 21.4 10000 167.5 10250 453.0 10500 935.7 10700 1500.4 10900 3118.1 11000 9283.3 11250 Infeasible 400 12250 12500 12750 13000 13250 13500 13750 6100. 7540. 9123. 11026. 13170. 24287. Infeasible 14000 19450. 14500 24434. 15000 29742. 15250 32079. 15500 34837. 15750 37966. 16000 64591. 16250 Infeasible 85 The intersection of the curve obtained from the results of suboptimization of a finite thermal capacity with the non-integrated operation line at F^ confirms the suboptimization process at this point. The rapid increase in thermal energy production and hence ther-mal fuel consumption when the finite capacity subsystem is operated at an annual firm energy output close to its maximum capability suggests that operation in this region is economically unattractive. Figure 8.2.1 also indicates that the optimal integrated operation of the finite thermal capacity subsystem results in the same output-cost performance as the infinite thermal capacity subsystem except at annual firm energy outputs close to the maximum attainable by the finite thermal capacity subsystem. Thus a close approximation to the finite capacity curves for a l l thermal capacity limits can be generated by calibrating the infinite thermal capa-city curve (Fig. 8.1.1) using Equation 8.2.2 to locate the vertical por-tion of the output-cost curve for each thermal capacity value. The results of calibrating a performance cost curve in this fashion are shown in Figure 8.2.2. The approximation described above, which only affects interpre-tation of the calibrated curve at annual energy outputs close to the limit for a given thermal capacity, can be allowed for by modifying Equations 8.2.1 and 8.2.2 as follows: T | = 12T - A 8.2.3 i F = F . + T 8.2.4 T I T where A represents a small fraction of the thermal capacity. 87 The effect of this modification is to confine interpretation of the infinite thermal capacity curve to the region where i t coincides exactly with the finite capacity curve. This is demonstrated in Figure 8.2.3. The results obtained for the Portage Mountain subsystem in-dicate that a value of A in the order of 600 mkwh will achieve the de-sired effect. This f u l l adjustment may not be appropriate when substan-t i a l reserve thermal capacity exists in the generating system. The calibration of air infinite thermal capacity curve as opposed to repeated suboptimization over a range of thermal capacity values results in a very significant computational reduction when produc-ing a set of finite thermal capacity output-cost curves. These curves are required in the co-ordination of the overyear storage subsystems as described in the following chapter. Present Worth Expec ted Cost mkwh ( the rma l energy ) 50 W I o > H td w a H z H W Z 3 S o W D M K •-3 •< O c o > f f o o M > H o z "3 M a W 03 to 88 CHAPTER 9 CO-ORDINATION OF THE LONG-TERM OPERATION OF THE HYDROTHERMAL SUBSYSTEMS AND SUBGROUPS The principal problem considered in this study is the long-term operation of a mixed generating system. The adoption of the hierarchical approach to the analysis of the. system lead fir s t to the decomposition of the system into subsystems, then to a demonstration of suboptimization for one of the subsystems, and finally in this chapter to the reco-ordination of the subsystem operation to achieve efficient integrated long-term operation of the system as a whole. 9.1 CO-ORDINATION OF OVERYEAR STORAGE SUBSYSTEMS The objective which will be imposed upon the overyear storage subgroup is to meet a prescribed firm energy output at minimum present worth expected cost subject to the constraint of the thermal generating capacity (again defined as maximum generating capability in any month and hence in any year) allocated to the subgroup. An example will be taken from the B.C. Hydro generating system which will shortly include tvjo major overyear storage projects: the Portage Mountain project and the Mica project. A similar analysis to that performed for the Portage Mountain project but using Mica project specifications, inflow data, etc. would provide similar results to those described in Chapter 8 for the Mica project and hence the output-cost curves 89 9 0 for the two projects is shown in Figure 9 . 1 . 1 . If F is again used to denote an annual firm energy output and C and T are the present worth expected cost and -annual thermal generating capacity requirement associated with F, and i f suffixes P, M, and O are used to denote the Portage Mountain, Mica project, and the subgroup formed by these two projects, then the co-ordination problem may be stated-Minimize CQ = C P + C M 9 . 1 . 1 Subject to the constraints that the sum of the subsystem firm energy outputs must equal a specified subgroup firm energy output. F = F + F 9 1 2 0 P M * * and that the thermal capacity allocation to the subgroup must be large enough to meet the requirements of the subsystems: TQ = Tp + TM 9 . 1 . 3 This problem can be solved graphically using the information con-tained in Figure 9 . 1 . 1 by plotting C_ against either F„ or Fw while meeting O P M the combined firm energy output prescribed by Equation 9 . 1 . 2 . Thus for a single value of FQ this will result in the graph shown in Figure 9 . 1 . 2 . * * If Fp and F^ are used to denote the contributions of firm energy * from each of the subsystems for minimum cost at subgroup output FQ and T and T^ the associated thermal capacity requirements then the thermal generating capacity required by the subgroup is given by: * * * T = T + T 9 1 3 O X P AM ' Portage Mountain FIGURE 9.1.1 ANNUAL FIRM ENERGY OUTPUT — COST CURVES FOR PORTAGE MOUNTAIN AND MICA SUBSYSTEMS 92 FIGURE 9.1.2 MINIMUM COST ALLOCATION OF FIRM ENERGY OUTPUT FOR PORTAGE MOUNTAIN AND MICA SUBSYSTEMS 93 But Equation 8.2.2 indicates that: and T = - F 9.1.4 P P P TM = FM'FM, 9'1-5 Therefore T0 = (FP " \> + (FM " FM l} " Fo " \ + FM l> 9 a'6 provided that F > (F„ + F.. ) otherwise T will be zero, o P^ 0 Equation 9.1.6 indicates that the thermal capacity allocation (as defined here) to the subgroup depends only upon the firm energy output of the subgroup and the maximum firm energy outputs of the overyear storage hydro projects when operated independently of any thermal plant. By repeating the subgroup analysis over a range of values for F^ and plotting the resulting minimum present worth expected cost values a sim-ilar curve to that obtained for the individual subsystems is obtained as shown in Figure 9.1.3. This analysis will produce a series of unique load alloca-* * tions (F , F ) over the combined output range: (F„ + F ) < F < (F + F.. ) Pl Ml - ° - P2 M2 which are shown graphically in Figure 9.1.4. Above and below this range of combined firm energy outputs the subsystem output-cost curves have identical slopes and the allocations are not uniquely defined. FD +FM Fo +F, M i FIGURE 9.1.3 ANNUAL FIRM ENERGY OUTPUT — COST CURVE FOR THE OVERYEAR STORAGE SUBGROUP 95 MINIMUM COST ANNUAL FIRM ENERGY ALLOCATIONS OVER A RANGE OF OVERYEAR STORAGE SUBGROUP OUTPUTS 96 The minimum cost firm energy output allocation curve (Figure 9.1.4) provides a basis for determining joint operating policies for the two overyear storage subsystems and hence the two overyear storage hydro projects when the combined firm energy output requirement is known. A procedure for establishing this output by co-ordinat-ion of the overyear storage subgroup and the rest of the system is discussed in the following section. 9.2 CO-ORDINATION OF OVERYEAR STORAGE AND SEASONAL STORAGE SUBGROUPS The optimization of the integrated operation of a number of single purpose seasonal storage hydro projects and thermal plants has been demonstrated in a number of publications [9] and suggests that the subopti-mization of the seasonal storage subgroup is entirely feasible. It is reason-able to assume that the optimum performance of this subgroup when measured in terms of minimum present worth expected thermal costs as before will be of the same form as was obtained for the overyear storage subgroup in Figure 9.1.3. The seasonal storage subgroup will be able to meet some firm energy output F without consumption of any thermal energy. Above this output Yl the cost of integrated operation will rise smoothly as shown in Figure 9.2.1 until some output F is reached. Any further increase in firm -energy output 2 beyond F will be obtained at an incremental cost eauivalent to independent 2 operation of the thermal plant. The procedure for allocating the total system firm energy demand F to the overyear storage and seasonal storage subgroups is similar to that performed for the overyear storage subsystems in the preceeding section. In this case the information contained in the graphs.shown in Figures 9.1.3 and (Postulated) 9 8 9 . 2 . 1 is used and the minimum cost firm energy output allocation performed at only one combined output F . Using suffixes S, 0 and Y to identify quantities associated with the overall system, overyear storage, and seasonal storage subsys-tems the subgroup co-ordination problem can be written! Minimize C = C + C 9 . 2 . 1 subject to* F = F + F 9 2 2 S * 0 Y Ts = TQ + TY 9 . 2 . 3 * If F Q and F ^ denote the minimum cost firm energy output allocations to the two subgroups then from Equation 8 . 2 . 2 : T S = ( F o " F 0 } + ( F Y " F Y } 9 * 2 * 4 = F - (P + F ) 9 . 2 . 5 S ° 1 Y l provided that F > (F + F ), otherwise T will be zero. 1 1 S This again suggests that the thermal capacity allocation as defined in this analysis is dependent only upon an overall firm energy output and the maximum firm energy outputs of the subgroups (and hence subsystems) when operated without thermal generation. 99 9.3 RECOVERY OF THE OVERYEAR STORAGE HYDRO PROJECT OPERATING POLICIES When the minimum cost firm energy allocations have been deter-mined for the major subgroups then the appropriate allocations to the in-dividual subsystems can be obtained, in the case of the overyear storage subsystems, from Figure 9.1.3. The operating policy for an individual over-year storage hydro project can then be determined by repeating the suboptimi-zation procedure at the appropriate subsystem annual firm energy output. The resulting optimal reservoir state trajectories (see Section 4.2) will provide details of monthly discharge decisions for a l l discrete annual inflow volumes and will constitute the operating policy for that hydro project within the overall generating system. An alternative and more refined approach, particularly in the case of the overyear storage projects, utilizes the state values obtained from suboptimization at the ^minimum cost firm energy output allocations in conjunction with a Dynamic Programming analysis of the combined operation of the two overyear storage projects over a one year period of operation. The state values serve as terminal state values at the end of the one year combined operation analysis and thereby resolve the unbounded nature of the problem beyond the current year of operation (see Section 3.3), The Dynamic Programming analysis of combined operation of the two hydro projects involves two state variables representing the storage volume in each reservoir at the beginning of a monthly time stage. The analytical refinemenet of this two state variable approach is attractive. However, the computational load increases exponentially with the number of state variables and consequently is impractical for the co-ordination of more than two projects even over the limited number of stages involved here. 100 9.4 CONCLUDING SUMMARY In this thesis the problem of determining an efficient policy for the long-term integrated operation of overyear storage hydro projects in a mixed hydroelectric and thermal electric generating system is solved by a system decomposition approach. An analytical method was fi r s t developed for determining the optimum integrated operation of a subsystem which includes an overyear storage hydro project and a thermal plant. The analysis combined a de-tailed Dynamic Programming optimization of the operation of this subsystem over a one year period together with a Markov decision process analysis of the long-term annual operation. Uncertainty regarding future annual in-flows and the unbounded time horizon of the decision process were both re-solved in this analysis. An application of the analysis was demonstrated for a subsystem which included the Portage Mountain hydro project. The results provide long-term integrated operating policy at a specified annual firm energy out-put and are directly applicable in the case of a generating system which in-cludes a single overyear storage hydro project, as is the present case in the B.C. Hydro system. The subsystem analysis was found to be computationally efficient and permitted repetition of the analysis of the Portage Mountain subsystem over a range of annual firm energy outputs and thermal capacity allocations. The results obtained were plotted to produce an energy output/cost curve for the subsystem. This curve was used as the basis of a minimum cost co-101 ordination procedure for the long-term operation of two overyear storage subsystems. A similar procedure, which employed the results obtained from the co-ordination of the two overyear storage subsystems, was then demon-strated for the co-ordination of these subsystems with the remainder of the generating system and provided an overall system operating policy. The individual overyear storage hydro project operating decisions associated with an overall system operating policy are retraced by repeating the subsystem analysis at the subsystem annual firm energy output prescribed by the policy. This analysis also produces a set of reservoir state values, appropriate to the overall policy, which may be utilized as terminal state values in a more detailed short-term analysis incorporating, for example, the direct interaction between two overyear storage hydro projects. In a planning context the analysis and co-ordination procedure described here could be performed for prospective overyear storage projects. The results would indicate the value of the proposed addition to the long-term integrated operation of the generating system and could also be used in the determination of the most advantageous ratio of thermal generating capacity to overyear storage capacity in the future development of the system. R E F E R E N C E S [1] Buras, Nathan. Scientific Allocation of Water Resources. New York: American Elsevier Environmental Science Series, 1972. [2] Silver, R., M. Okun and S. O. Russell. "Dynamic Programming in a Hydro-electric System," in Proceedings of the International Symposium on Model-ling Techniques in Water Resources Systems (Ed. A. K. Biswas), Ottawa, May, 1972, Vol. 2, p. 423. [3] Hall, W. A., A. J . Askew, and W. W. Yeh. "Use of the Critical Period in Reservoir Analysis," Water Resources Research, Vol. 5 (December, 1969), pp. 1205-1215. [4] Hufschmidt, M. M., and M. B. Fiering. Simulation Techniques for Design of Water Resource Systems. Cambridge, Mass.: Harvard Univer-sity Press, 1966. [5] Gablinger, M., and D. P. Loucks. "Markov Models for Flow Regulation," Proceedings, ASCE, Hydraulics Division, Vol. 96 (January, 1970), pp. 165-181. [6] Hall, W. A., W. S. Butcher, and A. Esogbue. "Optimization of the Oper-ation of a Multiple Purpose Reservoir by Dynamic Programming," Water Resources Research, Vol. 4 (June, 1968), pp. 471-477. [7] Hall, W. A., G. W. Tauxe and W. W. Yeh. "An Alternative Procedure for the Optimization of Operation for Planning With Multiple River Program Systems," Water Resources Research, Vol. 5 (December, 1969), pp. 1367-1372. [8] Parikh, S. C. "Linear Dynamic Decomposition Programming of Optimal Long-Range Operation of a Multiple Multi-Purpose Reservoir System," Proceed-ings, 4th International Conference on Operations Research, Boston, September, 1966. [9] Hall, w. A., and J . A. Dracup. Water Resources Systems Engineering. New York: McGraw-Hill CO. Inc., 1970. [10] Windsor, J.S., and Ven Te Chow. "Multi-Reservoir Optimization Model," Proceedings, ASCE, Hydraulics Division, Vol. 98 (October, 1972), pp. 1827-1844. 102 103 [11] Nemhauser, G. L. Introduction to Dynamic Programming. New York: John Wiley, 1970. [12] Hammer, P. C. (Ed.). Advances in Mathematical Systems Theory. Univer-sity Park, Penn.: Pennyslvania State University Press, 1969. [13] Mesarovic, M. D., D. Macko and Y. Takahara. Theory of Hierarchical and Multi-Level Systems. New York: Academic Press, 1970. [14] Arventidas, N. V., and J. Rosing. "Composite Representation of a Multi-Reservoir Hydro Power System," in Transactions, IEEE, Power Apparatus and Systems Division, Vol. 58 (January, 1970). [15] Wagner, H. M. Principles of Operation Research. Englewood C l i f f s , N.J.: Prentice-Hall, 1969. [16] Conard, J . W. An Introduction to the Theory of Interest. Berkeley: University of California Press, 1959. [17] Parzen, E. Stochastic Processes. San Francisco: Holden Day, 1962. [18] Howard, R. A. Dynamic Programming and Markov Processes. Cambridge, Mass.: M.I.T. Press, 1960. [19] Chow, Ven Te. Handbook of Applied Hydrology. New York: McGraw-Hill . Co. Inc., 1964. [20] Elderton, Sir W. P., and N. L. Johnson. Systems of Frequency Curves. London: Cambridge University Press, 1969. [21] Private correspondence, B.C. Hydro and Power Authority, Vancouver, 1971. [22] Young, G. K. "Finding Reservoir Operating Rules," Proceedings, ASCE3 Hydrology Division, Vol. 93 (November, 1967), pp. 297-321. [23] Smillie, K. W. An Introduction to Regression and Correlation. Toronto: Ryerson Press, 1966. [24] Larson, Robert E. State Increment Dynamic Programming. New York: American Elsevier Publishing Co. Inc., 1968. A P P E N D I C E S APPEND LX 1 COMPUTER PROGRAM LISTING 1 C MAIN PROGRAM • 2 C 3 DIMENSION V (20) ,P (20,20) ,0 (20) 4 DOUBLE PRECISION BB (20) ,AA (20,20) ,DET 5 DIMENSION TEMP (20,20) ,PP (20,20) 6 DIMENSION PHOBIN (9) ,X (8) , Y (8) , A (1 3) , B (13) ,EMONTH (13) 7 DIMENSION WAV1H(9) 8 COMMON EUPPER,ELOWER,NSTATE,STINC,PROBIN,X,Y,A,E,EMONTH 9 1 ,A1IEDEM, SYS AUG, KTRIG,WAVIN, DISC, I PEAS 10 c 11 c 12 DISC=0.926 13 KK=0 14 ITRLIM=10 15 NSTATE=20 16 EUPPER=2200. 17 ELOWER=2155. 18 STINC=(EUPPER-ELOWER)/ (NSTATE-1) 19 c SPECIFY ANNUAL FIRM ENERGY OUTPUT ANEDEM 20 ANEDEM=12000. 21 DEKLIM=12600. 22 DEMINC=200. 23 c SPECIFY MONTHLY THERMAL GENERATING CAPACITY 24 SYSAUG=300. 25 c MTRIG=1 FOR FULL DP PRINTOUT, =0 IF NOT 26 MTRIG=0 27 READ (5,9) (WAVIN (I) , PROBIN' (I) ,1=1,9) 28 9 FORMAT (2F10.3) 29 READ (5,10) (X (I) ,Y (I) ,1=1,8) 30 10 FORMAT (2F10.3) 31 READ (5,12) (A (ISTAGE),B (ISTAGE) ,ISTAGE=1,13) 32 12 FORMAT (2F10.5) 33 READ (5, 13) (EMONTH (ISTAGE) ,ISTAGE=1, 13) 34 13 FORMAT (13F5.3) 35 c INITIALISATION OF RESERVOIR STATE VALUES V (J) 36 V (1) =1.E8 37 DO 5 J=2,NSTATE 38 V (J) =V (J-1)-3300. 39 IF (J.EQ.10)V(J)=66400. 40 5 CONTINUE 41 19 ITER=0 42 c 43 c 44 c 45 c ITERATIVE CYCLE COMMENCES 46 20 ITER=ITER•1 47 IF(ITER.LE.ITRLIM)GO TO 22 48 KRITE(6,21) ITRLIM 49 21 FORMAT('CONVERGENCE NOT ACHIEVED IN ',12,' CYCLES') 50 GO TO 204 51 c 52 c POLICY IMPROVEMENT BY DYNAMIC PROG. SBRTN 53 22 CALL DP 1 (V,P,Q) 54 c 55 c VALUE DETERMINATION 106 107 56 DO 15 1=1,NSTATE 57 U»(I)=-Q(I) 58 DO 16 J=1,NSTATE 59 AA (I,J) = DISC* T (I, J) 60 IF(I.EQ.J)AA(I,J)=AA(L,J)-1.0 61 16 CONTINUE 62 15 CONTINUE 63 CALL DSOLTN (AA , F.B, NSTATE, NSTATE, DET) 61 IF(DET)201,201,203 65 201 WRITE (6,202) 66 202 FORMAT('INVERSION FAILED - PROGRAH TERMINATED') 67 GO TO 204 68 C 69 c TEST FOR CONVERGENCE OF V(J) 70 203 N=0 71 WRITE (6,25) 72 25 FORMAT (' V (I) BB(I) Q(I)1) 73 DO 17 1=1,NSTATE 74 WRITE (6, 18) V (I) , BB (I) ,Q (I) 75 18 FORMAT (2X,3F12.1) 76 IF(BD (I) .GT.1.E5)GO TO 130 77 IF (DABS (BB (I)-V (I)) .LT. 1.) N=N*1 78 GO TO 132 79 130 IF (ITER. GE. 2) N=HM 80 132 V(I)=BB(I) 81 17 CONTINUE 82 IF (N.LT.NSTATE)GO TO 20 83 122 IF (MTRIG.11E.2) GO TO 100 84 WRITE(6,120) 85 120 FORMAT(///20X,'OPTIMAL DP TRAJECTORIES'//) 86 WRITE (6,121)ANEDEM,SYSAUG €7 121 FORMAT ('ANNUAL ENERGY DEMAND =•,F10. 1,1X, 'MKWH•,5X, 88 1'SYSAUG =',F10.1,» MKWH'///) 89 GO TO 20 90 C 91 C 92 C 93 100 WRITE (6,105) 94 105 95 IF (IFEAS.EQ.1)GO TO 106 96 WRITE(6,107) 97 107 FORMAT {/• AT LEAST ONE INFEASIBLE TRAJECTORY'/) 98 106 WRITE (6, 104) ( (P(I,J) ,J=1, NSTATE) ,1=1, NSTATE) 99 104 FORMAT(20F6.3) 100 WRITE(6,103) (V (I) ,Q (I) ,1=1,NSTATE) 101 103 FORMAT (2F15.1) 102 C 103 c CALCULATION OP STEADY STATE PROBY. MATRIX 104 MM=0 105 DO 605 1=1,NSTATE 106 DO 606 J=1,NSTATE 107 PP(I,J)=P(I,J) 108 606 CONTINUE 109 605 CONTINUE 110 600 MM=MH+1 111 CALL MULT (P,PP,TEMP,NSTATE,NSTATE) 112 616 DO 601 1=1,NSTATE 113 DO 602 J=1,NSTATE 114 IF (ABS(TEMP(I,J)-PP (I,J)).GT.0.001)GO TO 610 115 602 CONTINUE 108 116 601 CONTINUE 117 GO TO 620 118 610 DO 603 1=1,NSTATE 119 DO 604 J=1,NSTATE 120 PP (I,J)=TEMP (I,J) 121 604 CONTINUE 122 603 CONTINUE 123 IF(MM.LT.20)GO TO 600 124 620 WRITE (6,613)KM 125 613 FORMAT ('STEADY STATE PROBABILITIES. ',13,' TR At.':. ITIONS •) 126 615 WRITE (6,611) ((TEMP (I,J) ,0=1, 11 S T A T E ) ,1=1, NST AT E) 127 611 FORMAT (20F6.3) 128 PWEC=0.0 129 DO 630 J=1,NSTATE 130 PWEC=PWEC*V (J)*TEMP (I,J) 131 630 CONTINUE 132 WRITE (6,632)DISC 133 632 FORMAT (//• DISCOUNT FACTOR *,P10.3//) 134 WRITE (6,631)SYSAUG,ANEDEM,PWEC 135 631 FORMAT(//• MAX. AUGMENTING ENERGY AVAIL. IN ANY MONTR *, 136 1F10.1,» MKWH'//' ANNUAL FIRM ENERGY EEM.',F10.1, 137 2' MKWH'//• PRESENT WORTH EXPECTED COST OVER', 138 3» RECURRENT STATES',F10.1,' COLLARS'//) 139 ANEDEM=ANEDEb+DEMINC 140 IF (IFEAS.F.Q.O) GO TO 204 141 IF(PWEC.GT.1.E8JG0 TO 204 142 IF(ANEDEM.LT.DEMLIM+20.)GO TO 19 143 C 144 204 STOP 145 C 146 END 1> 1 •> C r » DYNAMIC PROGRAMMING SUBROUTINE FOR POLICY IMPROVEMENT 3 L . SUBROUTINE DPI (V,P,Q) 4 DIMENSION IFOL-(13,20) ,GOPT (13,20) 5 DIMEliSIOS F (20) ,PROHIN (9) , JEND (20) 6 DIMENSION COST (20) ,E (20) ,D (20) 7 DIMENSION FOPT (13,20) ,DOPT (13,20),ELX (20) ,Q1 (4) 8 DIMENSION V (20) ,P (20,20) ,Q (20) 9 DIMENSION X (8) ,Y (0),A (13),B (13),EMONTH (13) 10 DIMENSION PELX{13,20),PD(13,20),PG(13,20) 11 DIMENSION WAVIN (9) 12 COMMON EUPPER,ET.OWER,NSTATE,STIHC,PROBIN,X,Y, A, B, EMONTH 13 1,ANEDEM,SYSAUG,MTRIG,WAVIN,DISC,IFEAS 14 c 15 c 16 NDECIS=20 17 a=i 18 IFEAS=1 19 c INITIALISE P ( ) , Q() TO ZERO; E () TO STATE ELEVATIONS; 20 DO 306 1=1,NSTATE 21 DO 307 J=1,NSTATE 22 P(I,J)=0.0 23 307 CONTINUE 24 Q(I)=0.0 25 E(I)=(FLOAT(I)-1.)*STINC+ELOHER 26 306 CONTINUE 27 c 28 c 29 c DP ANALYSIS AND TRAJECTORY RECOVERY COMMENCES 30 DO 301 INFLO=1,9 31 c 32 c DP ANALYSIS FOR A PARTICULAR INFLOW COMMENCES 33 DO 320 1=1,NSTATE 34 FOPT(1,I)=DISC*V(I) 35 320 CONTINUE 36 DO 303 ISTAGE=2,13 37 IF (MTRIG.NE.2)GO TO 931 38 IF (IHFL0.LT.7)GO TO 931 39 IF (ISTAGE.NE.2)GO TO 931 40 WRITE{6,950) (E(I) , FOPT(ISTAGE-1,1),1=1,NSTATE) 41 950 FORMAT (2F12.1) 42 931 DO 304 ISTATE=1,NSTATE 43 FMIN=1.E12 44 DO 305 IDECIS=1,NDECIS 45 ELINIT= (FLOAT(ISTATE)- 1.)*STINC + ELOWER 46 IELFOL=ISTATE*IDECIS-4 47 IF(IELFOL.GT.NSTATE)GO TO 305 48 IF(IELFOL.LT.1)GO TO 305 49 ELFOL=FLOAT(IELFOL-1)*STINC+ELOWER 50 c 51 c REFER TO SUBPROGRAM SIM1 FOR SIMULATED OPERATION 52 CALL SIM1(ISTAGE,ELI NIT,ELFOL,INFLO,G,DISCHF) 53 c 54 c IF STATE CHANGE GIVES NEC. DISCH. TEST NO LARGER 55 IF(DISCHF.LT.0.0)GO TO 315 109 110 56 FF=G•FOPT (ISTAGE- 1 , IELFOL) 57 IF (KTRIG.NE.2)GO TO 930 58 IF (INF1.0. LT. 7) GO TO 930 59 IF (ISTAGE. HE. 2) GO TO 930. 60 WHITE (6,951) EL I NIT, ELFOL,DISCIIF, FOPT (ISTAGE-1 , IELFOL) 61 951 FORMAT (UF12. 1) 62 930 IF (FF.LT. FMIN) DECOPT=DISCHF 63 IF (FF. LT. PMIN) IELOPT=IELFOL 64 IF (FF.LT.FHIN)GBEST=G 65 IF (FF.LT.FMIN)FKIN=FF 66 305 CONTINUE 67 315 FOPT (ISTAGE,ISTATE)=FMIN 68 GOPT (ISTAGE,ISTATE)=GBEST 69 IFOL (ISTAGE,ISTATE)=IELOPT 70 DOPT (ISTAGE,ISTATE)=DECOPT 71 304 CONTINUE 72 303 CONTINUE 73 C 74 IF (MTRIG.NE.2)GO TO 311 75 309 DO 960 1=1,NSTATE 76 WRITE(6,961) (FOPT (14-11,1) ,11=1,12) 77 WRITE (6,961) (GOPT (14-11,1) ,11=1, 12) 78 WRITE(6,961) (DOPT (14-11,1) ,11=1,12) 79 961 FORMAT (12F10.1) 80 WRITE (6,962) (IFOL (14-11,1) ,11=1,12) 81 962 FORMAT (12110/) 82 960 CONTINUE 83 C 84 C TRAJECTORY RECOVERY 85 C 86 308 IF(MTRIG.NE.2)GO TO 311 87 WRITE (6,312)INFLO 88 312 FORMAT (1H1,20X,* DISCRETE TOTAL ANNUAL INFLOW',14//) 89 311 ISTAGE=13 90 C INITIALISE ELX() AND COST ACCUMULATOR COST () ARRAYS 91 DO 505 I=1,NSTATE 92 PELX (ISTAGE,I)=E (I) 93 COST (I)=0.0 94 505 CONTINUE 95 C RECOVERY CYCLE BEGINS 96 C IDENTIFY FOLLOW IMG STATES FOR ALL STATES 97 DO 503 ISTART=1,NSTATE 98 ISTATE=ISTART 99 ISTAGE=13 100 501 COST (ISTART)=COST(ISTA RT)+GOPT(ISTAGE,ISTATE) 101 PD (ISTAGE,ISTART)=DOPT (ISTAGE,ISTATE) 102 PG (ISTAGE,ISTART)=GOPT (ISTAGE,1STATZ) 103 PELX (ISTAGE,ISTART)=E (ISTATE) 104 ISTATE=IFOL(ISTAGE,ISTATE) 105 ISTAGE=ISTAGE-1 106 IF (ISTAGE.GT.1)GO TO 501 107 JEND (ISTART)=ISTATE 108 503 CONTINUE 109 IF(MTRIG.NE.2)GO TO 510 110 DO 511 ISTATE=1,NSTATE 111 WRITE(6,512) (PELX(14-11,ISTATE),11=1, 13) 112 512 FORMAT (13F10. 1) 113 WRITE (6,513) (PD (14-11, ISTATE) ,11=1, 13) 114 513 FORMAT(13F10.1) 115 WRITE (6,514) (PG (14-11,ISTATE) ,11=1, 13) I l l 116 514 FORMAT(13F10.1//) 117 511 CONTINUE 118 C TRAJECTORIES, ENDSTATHS, AND COT.TS RECOVEREC FOR ALL 119 c 120 c ACCUMULATE P(I,J) AND Q(I) 121 510 DO 506 1=1,NSTATE 122 P(J,JE!lD(I))=P(I,JEND(I))+PROniN (IN FLO) 123 Q (I) ~Q (I) + COST (I) *PROIIIH (INFLO) 124 506 CO NTINUE 125 301 CONTINUE 126 c 127 c 128 940 WRITE(6,905) ( (P(1,0) ,0=1,NSTATE) ,1=1,NSTATE) 129 905 FORMAT (20F6.3) 130 GO TO 399 131 908 WRITE (6,906) (Q(I) ,1=1, NSTATE) 132 906 FORMAT (F10.1) 133 WRITE(6,907) (V (I) ,1=1,NSTATE) 134 907 FORMAT (F10.1) 135 GO TO 399 136 350 WRITE (6,351)ISTAGE,INFLO,ELFOL 137 351 FORMAT ('ELEVATION AT END OF STAGE EXCEEDS LIMITS 138 1. PROGRAM STOPPED',216,F10.1) 139 c . ; — 140 399 RETURN 141 c 142 END 1 2 C C SIMULATION SUBROUTINE FOR MONTHLY DISCHARGE DECISION 3 4 C SUBROUTINE SIK1 (ISTAGE,ELINIT,ELFOL,INFLO,G,DISCHF) 5 DIMENSION X ( 8 ) , Y ( 8 ) , Q 2 (6) 6 DIMENSION A (13) ,B (13) ,EMONTH (13) 7 DIMENSION PROBIN (9) ,H AVIN(9) 8 COMMON EUPPER,SLOWER,NSTATE,STINC,PROBIN,X,Y,A,E,EMONTH 9 1,AHEDEM,SYSAUG,MTRIG,HAVIN,DISC, IFEAS 10 C 11 C 12 M=2 13 ETAIL=1649.0 14 EFFY=0.9 15 C 16 C COMPUTE DISCHARGE VOL. AND FLOW, MEAN HEAD 17 VIN= (A (ISTAGE)+ B (ISTAGE)*WAVIN(IHFLO))*1.E9 18 VOL=SAINT (8,X,Y,ELIHIT,M,Q2) 19 VOLFOL=SAIIJT (8 , X, Y, ELFOL , H, Q2) 20 DISCHV=VIN+(VOL-VOLFOL)*43560.*1.E6 21 IF (DISCHV.LT.O.O)DISCHF=-1.0 22 IF (DISCHV.LT.0.0)GO TO 801 23 DISCHF=DISCHV/ (730.*3600.) 24 HMEAN= ( (ELINIT + ELFOL)/2.0)-ETAIL 25 DISHAX=73000.0- ( (HMEAN+ETAIL-2170.) *7000./30.) 26 DISTEM=7300 0.-(2170.-(HMEAN+ETAIL))*3000./30.0 27 IF (HMEAN+ETAIL.LT.2170.)DISMAX=DISTEM 28 IF(ELFOL.GT.2199.0)DISMAX=1.E6 29 IF(DISCHF.GT.DISMAX)GO TO 801 30 C 31 C ENERGY GENERATED BY DISCHV IN MKWH 32 ENGY=DISCHV*HMEAN*EFFY*0.235*1.E-10 33 EAUG=EMONTH (ISTAGE) *ANEDEM-ENGY 34 IF(EAUG.LT.O.O)EAUG=0.0 35 G=EAUG 36 IF (EAUG.LE.SYSAUG)GO TO 800 37 801 G=1.E12 38 C 39 800 RETURN 40 C 41 END 112 APPENDIX 2 GENERATION OF ANNUAL INFLOW VOLUME FROM GENERATED MONTHLY INFLOWS Let: represent a recorded monthly inflow from month m, year i , where m = 1,. . .,12, i = 1,. . . , N represent a recorded annual inflow in year i , The simple regression model (Equation 7.3.1) used to generate an inflow volume in month y from an annual inflow volume X i s : m ym " bom + blm( X ) f o r m = 1" ' "1 2 where: b £ - -.(x.-x)(y -y ) x x xm m lm „ , -.2 E(x^-x) i Numerator: = Ex.y. - Ex.y - Exv + Exv . i im . x'*n . -'im . I i m . m Ex.y. - Nxy - Nxy + Nxy . a/xm m m ;m Ex.y. - Nxy . x xm m x 113 D e n o m i n a t o r : E ( x i - x ) 2 H e n c e : 2 -2 E x - N x l E b m E [ E x . y . - N x y ] . . „ . i i m .m _ m I ' l m ~ 2 -2 E x . - N x . i x N u m e r a t o r : E E x . y . - N x E y m x m T h e r e f o r e : E x . E y . - N x E E ^ y . A x J x m m x N - ' x m E x . x . - x E x . . X X . X E x ? - N x 2 X E b m l m T 2 4x. x x N x 2 E x 7 - N X = 1.0 A l s o : b = y - b , x o m m l m 115 Therefore: Eb = E[y - b, x] om m lm m m = x - xEb, lm = X - x 0.0 For any given annual inflow volume X the synthesized total X will be the sum of the generated monthly inflows y ; m = 1,. . .,12 x = Eym m m E(b + b. (X)) „ om lm m = Eb + XEb, om lm m m which, from the previous page, has been shown to be: = 0.0 + 1.0X = X A Hence X = X so that the generated annual total inflow will always equal the given value.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The long-term operation of overyear storage hydroelectric...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The long-term operation of overyear storage hydroelectric projects within a mixed hydrothermal generating… Caselton, William F. 1973
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | The long-term operation of overyear storage hydroelectric projects within a mixed hydrothermal generating system |
Creator |
Caselton, William F. |
Publisher | University of British Columbia |
Date Issued | 1973 |
Description | It is necessary to establish a new operating policy when a major hydro electric project is introduced into an existing generating system. When the project has a large reservoir, capable of storing substantial volumes of water over a number of years, then the operating policy cannot be confined to the current years operation - a long term operating policy is required. Analysis of the operation of the hydro project to determine the optimum long term policy becomes complex when the benefits of complementary interaction of the project with the rest of the generating system are considered. Further complication arises in the analysis from the uncertainty associated with reservoir inflows over this long time period. These complications prohibit the determination of an optimal long term operating policy by direct application of presently available optimisation techniques. In this thesis a system decomposition approach is presented for the determination of efficient long term operating policies for one or more overyear storage hydro projects in a mixed hydro and thermal electric generating system. The system is first decomposed into a set of subsystems which are designed to preserve the complimentary operation of individual hydro and thermal plants. The long term operation of a subsystem is optimised over a range of annual firm energy outputs using an analytical procedure which simultaneously applies optimisation techniques to both short term and long term operation. When only one overyear storage project exists in a generating system the results obtained from the analysis of its associated subsystem provide the required long term operating- policy. In the case where a number of overyear storage projects exist, or are proposed, it is necessary to first analyse the associated subsystems individually. The results obtained form the basis of a procedure for determining an efficient operating policy which, over the long term, coordinates the operation of the overyear storage projects and at the same time integrates their operation with thermal and short term storage plants within the generating system. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-25 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050537 |
URI | http://hdl.handle.net/2429/31828 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1973_A1 C38_4.pdf [ 5.22MB ]
- Metadata
- JSON: 831-1.0050537.json
- JSON-LD: 831-1.0050537-ld.json
- RDF/XML (Pretty): 831-1.0050537-rdf.xml
- RDF/JSON: 831-1.0050537-rdf.json
- Turtle: 831-1.0050537-turtle.txt
- N-Triples: 831-1.0050537-rdf-ntriples.txt
- Original Record: 831-1.0050537-source.json
- Full Text
- 831-1.0050537-fulltext.txt
- Citation
- 831-1.0050537.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050537/manifest