GENETIC ALGORITHMS FOR MULTI-OBJECTIVE OPTIMIZATION IN WATER QUALITY MANAGEMENT UNDER UNCERTAINTY by BRYAN ANTONY TOLSON B.Sc. (Env) The University of Guelph, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Civil Engineering Department) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 2000 Â© Bryan Antony Tolson, 2000 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g of t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the head of my department or by h i s or her r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d without my w r i t t e n p e r m i s s i o n . Department The U n i v e r s i t y of B r i t i s h Columbia Vancouver, Canada Date ABSTRACT This thesis demonstrates the combined usage of a number of novel approaches and techniques in the multi-objective management of water quality systems under uncertainty. The First-Order Reliability Method (FORM) is used to estimate the risk-based system performance indicators of reliability, vulnerability, and resilience that provide measures of the frequency, magnitude and duration of the failure of water resource systems, respectively. FORM accuracy and efficiency for performance indicator estimation is compared extensively with Monte Carlo Simulation (MCS). Genetic Algorithms (GAs) are demonstrated as a robust optimization technique by solving various multi-objective water quality management models that optimize the performance indicators and the total point source waste treatment cost. In addition, the Tradeoff Surface Representation (TSR) Algorithm is incorporated as a general multi-objective technique for accurate and efficient identification of convex tradeoff surfaces. The Willamette River Basin in Oregon, USA is utilized as the water quality management case study for the demonstration of all techniques. The performance indicators are estimated with respect to meeting dissolved oxygen (DO) standards and ambient DO is simulated using a QUAL2E water quality response model. Results show that FORM estimates of the performance indicators, while significantly less accurate than MCS estimates, seem to provide reasonable results when utilized within the multiobjective water quality management models. A comparison of FORM and MCS shows that while FORM is more efficient relative to MCS, the difference in efficiencies is significantly less than previously reported in the literature. The TSR Algorithm, in comparison with the commonly used Constraint Method for multi-objective tradeoff curve generation, is shown to produce a superior representation of the tradeoff curve. Furthermore, the TSR Algorithm is also shown to produce a maximum amount of auxiliary information regarding the bounds on the location of the tradeoff curve between tradeoff points. 11 TABLE OF CONTENTS ABSTRACT ii LIST OF TABLES vi LIST OF FIGURES viii ACKNOWLEDGEMENTS x 1 INTRODUCTION 1 2 LITERATURE REVIEW 5 2.1 Water Quality Management 5 2.2 Uncertainty in Environmental Response Modelling 6 2.3 Risk-Based System Performance Indicators in Water Resources 8 2.4 First-Order Reliability Method ( F O R M ) in Water Resources 10 2.5 Genetic Algorithms (GAs) in Water Resources 11 2.6 Multi-Objective Optimization 12 2.7 Recent Management Models for Water Resource Systems Under Uncertainty 16 3 METHODS OF ANALYSIS 3.1 Reliability Analysis 3.1.1 3.2 .' FORM 20 20 21 F O R M and Monte Carlo Simulation (MCS) Estimators of Reliability, Vulnerability and Resilience 23 3.2.1 Reliability 23 3.2.2 Vulnerability 23 3.2.3 Resilience 25 3.2.4 Accuracy and Efficiency of FORM and MCS Estimates of the Performance Indicators 28 3.3 GAs 31 3.3.1 Decision Variable Coding:....' 33 3.3.2 Initialization 33 3.3.3 Fitness Evaluation 35 3.3.4 Stopping Criteria 35 3.3.5 Selection 35 3.3.6 Crossover 36 3.3.7 Mutation 36 iii 3.3.8 3.4 4 Other Advanced GA Operators 37 Tradeoff Surface Representation (TSR) Algorithm 38 3.4.1 The General TSR Algorithm 39 3.4.2 Mathematical Formulation of the General m-Objective TSR Algorithm 45 3.4.3 Previous Similar Attempts at Efficient Tradeoff Curve Approximation 47 3.4.4 Summary Notes on the TSR Algorithm 50 W A T E R Q U A L I T Y M A N A G E M E N T IN T H E W I L L A M E T T E R I V E R B A S I N 4.1 52 Willamette River Basin and Water Quality System Response Model 4.1.1 4.2 52 Water Quality System Response Model 54 Estimation of Reliability, Vulnerability and Resilience 4.2.1 55 System Uncertainty 55 4.2.1.1 Natural Variability: Flow and Temperature Values 56 4.2.1.2 Parameter Variability: Sediment Oxygen Demand (SOD) Coefficient 60 4.2.1.3 Parameter Variability: Reaeration (Ka) Coefficient 62 4.2.1.4 Correlations Between Random Variables 66 4.2.2 Application of FORM to a Steady-State Response Model for Resilience Estimation 69 4.2.3 Location and Corresponding Critical Dissolved Oxygen (DO) Levels 71 4.2.4 Wasteload Levels 72 4.2.5 Estimation of the Performance Indicators 73 4.3 Multi-Objective Optimization of the Risk-Based Performance Indicators and Total Waste Treatment Cost 4.3.1 Waste Treatment Cost Data Development of Wastewater Treatment Plant (WWTP) Cost Data 75 4.3.1.2 Development of Pulp and Paper Mill (PPM) Cost Data 79 Computational Efficiency Improvements 4.3.2.1 4.3.3 4.3.3.2 4.3.4 82 Random Variable Reduction 83 Multi-Objective Water Quality Management Model Formulations 4.3.3.1 84 MCS as the Performance Indicator Estimation Technique in the Water Quality Management Models 86 GA Parameters 87 Method of Evaluating the Water Quality Management Solutions R E S U L T S A N D DISCUSSION 5.1 75 4.3.1.1 4.3.2 5 74 Application of F O R M and M C S for Performance Indicator Estimation 88 92 92 5.1.1 Accuracy of the FORM Estimates of the Performance Indicators Relative to MCS 93 5.1.2 Complete Assessment of FORM versus MCS for Estimates of the Performance Indicators 97 5.1.3 Additional Investigations 105 iv 5.1.3.1 Impact of Reliability and Vulnerability Correlation Matrix Modification 105 5.1.3.2 Effect of Parameter Auto-Correlation 106 5.1.3.3 FORM Non-Convergence and Design Point Interpretation 106 5.1.4 5.2 Summary 109 Application of the G A to Evaluate the Multi-Objective Water Quality Management Models 110 5.2.1 6 Preliminary Analyses for Specification of the Optimization Models 110 5.2.1.1 Random Variable Reduction 110 5.2.1.2 GA Sensitivity Analysis 113 5.2.2 Application of the TSR and Constraint Methods for Multi-Objective Optimization 115 5.2.3 FORM and MCS Comparison in the Cost-Reliability TSR Model 119 5.2.4 Tradeoff Curve Comparison Between Cost and Performance Indicators 123 CONCLUSIONS AND FUTURE WORK 129 NOTATION 136 BIBLIOGRAPHY 139 APPENDIX I. Stationarity Plots for the Data used in the Resilience Estimates 147 APPENDIX II. Waste Treatment Cost Data 149 APPENDIX III. Correlation Matrices used in the Water Quality Management Models : 151 V LIST OF TABLES Table 4.1. Information at utilized flow and temperature stations 57 Table 4.2. Annual seven-day moving average low flow statistics of flows and temperature used in reliability and vulnerability estimation 59 Table 4.3. 13-Day independent average statistics of flow and temperature data used for the resilience calculation 61 Table 4.4. Statistics for SOD (Tetra Tech 1995a) and Ka prediction errors used for all performance indicator estimation 63 Table 4.5. Original correlation coefficients calculated from the annual seven-day moving average low flow and temperature data used for reliability and vulnerability estimation 67 Table 4.6. Modified correlation coefficients used for the first stage reliability and vulnerability estimation 67 Table 4.7. Original correlation coefficients calculated from the 13-day independent average flow and temperature data used for resilience calculation Table 4.8. Modified correlation coefficients used for first stage resilience estimation 68 68 Table 4.9. Information about the failure states used for vulnerability estimation (ODEQ 1995) 72 Table 4.10. Summary of managed point source information 76 Table 4.11. Summary of economic factors required to update WWTP cost estimates to 1999 values 78 Table 4.12. Summary of economic factors required to update PPM cost estimates to 1999 average values 80 Table 5.1. FORM and MCS relative efficiency determination for reliability 102 Table 5.2. FORM and MCS relative efficiency determination for resilience 104 Table 5.3. Summary of the FORM and 5000-realization MCS estimates of the performance indicators for the reliability-, vulnerability- and cost-resilience tradeoff curve solutions Table II-1. Waste treatment cost data for the WWTPs and the PPMs 128 149 Table III-1. Original correlation coefficients calculated from the flow and temperature data used for reliability and vulnerability calculation in the water quality management models 151 vi Table III-2. Original correlation coefficients calculated from the flow and temperature data used for resilience calculation in the water quality management models 151 LIST OF FIGURES Figure 3.1. FORM approximation of the failure surface in standard normal space 22 Figure 3.2. Cumulative distribution function (CDF) of load (L) and schematic representation of multiple failure states for vulnerability estimation 25 Figure 3.3. Simple GA flow chart 34 Figure 3.4. Single-point crossover 36 Figure 3.5. Graphical representation of the TSR Algorithm; (a) step i (b) step ii (c) step iii 40 Figure 3.6. Determination of the tradeoff curve bounds using available auxiliary information in (a) step (iv) of the TSR Algorithm and (b) with the Constraint Method 42 Figure 3.7. Development of the mean estimate of the approximate tradeoff curve 44 Figure 3.8. Graphical presentation of a two-objective TSR Algorithm mathematical formulation based on Equation (32) 47 Figure 4.1. Willamette River Basin, Oregon, USA 53 Figure 4.2. General optimization program flow chart 90 Figure 5.1. FORM and MCS reliability estimates at RK 0 for a DO standard of (a) 5 mg/L and (b) 6 mg/L 94 Figure 5.2. FORM and MCS vulnerability estimates at RK 0 for two DO standards 95 Figure 5.3. FORM and MCS resilience estimates at RK 0 for a DO standard of (a) 5 mg/L and (b) 6 mg/L 96 Figure 5.4. Approximate FORM and MCS reliability operating curves for an actual reliability of (a) 0.5 and (b) 0.9 (or 0.1) 99 Figure 5.5. FORM and MCS reliability operating curves to assess the relative efficiencies of FORM versus MCS when the actual reliability is equal to approximately 0.5 102 Figure 5.6. Impact of modified correlations on FORM and MCS reliability estimates for a standard of 5 mg/L 106 Figure 5.7. Resilience at a DO standard of 5 mg/L with and without parameter auto-correlations 107 Figure 5.8. Normalized design points for time step 1 in resilience estimate at a DO standard of 5 mg/L with parameter auto-correlations for select random variables 108 Figure 5.9. Reliability curves for various stages of the random variable reduction process 112 viii Figure 5.10. Sensitivity analysis results for each of the best micro-GA and regular GA parameter sets 114 Figure 5.11. Cost-reliability tradeoff curves and partial tradeoff boundaries for the TSR and Constraint Methods 117 Figure 5.12. TSR and Constraint Method comparison at approximately the same reliability levels 118 Figure 5.13. TSR Algorithm tradeoff curves using FORM and 40 MCS realizations for reliability estimation 120 Figure 5.14. TSR Algorithm tradeoff curves using FORM, 40 and 80 MCS realizations for reliability estimation 122 Figure 5.15. TSR cost-vulnerability tradeoff curves for the FORM and 5000-realization MCS estimates of vulnerability 124 Figure 5.16. TSR cost-resilience tradeoff curves for the FORM and 5000-realization MCS estimates of resilience Figure I-1. Stationarity plot for the (a) mean and (b) standard deviation of the 13-day independent average flows over the low flow season Figure 1-2. 125 147 Stationarity plot for the (a) mean and (b) standard deviation of the 13-day independent average temperature over the low flow season 148 ix ACKNOWLEDGEMENTS First, I would like to thank my wife Kerry for giving me the extra strength I needed to keep working and for having the patience to put up with the long hours this work required. Many thanks to my supervisor, Dr. Barbara Lence, for all her patience, efforts and guidance regarding this thesis but also for her tremendous efforts, advice and discussions with respect to my future. Thanks to Dr. Ricardo Foschi for reading this work but most of all, for the collaboration that allowed this thesis to happen. Another important person without whom I would not have been able to accomplish as much as I did, is Dr. Hong Li who provided insights to the RELAN program, programming tips and countless modifications to the program in order for it to be utilized in this thesis. Thank you Dr. Holger Maier for providing some of the key ideas that started this thesis. Other important people that spent time answering questions, providing advice, sharing programs and data or assisted me in completing this work in any other way include James Bloom, Steve Schnurbusch and Mark Hamlin at the Oregon Department of Environmental Quality, Dr. David Carroll, Boris Fichot, Nicolas Lauzon, Melinna Schofield and Jo Miller. Finally, I would like to thank NSERC for providing me with the financial means necessary, in the form of a Postgraduate Scholarship, to pursue a Master's degree. 1 INTRODUCTION In surface water quality management, the classical problem is to determine the "best" set of waste treatment technologies for polluters, or the wasteload allocation, that achieves a specified water quality goal, or standard, and meets other societal goals such as minimal total pollution control costs. Effective environmental management requires a thorough understanding of the system being managed, suitable metrics to characterize all important aspects of system performance and a method for determining the management strategy that best meets the management objectives. Environmental response models are used to represent our state of knowledge about how the environment responds to anthropogenic impacts. They are mathematical representations that relate the rate of pollutant emission (i.e. mass or concentration of pollutant) at some location, and the resulting pattern of effects (i.e. the concentration of the pollutant) at various locations in the receiving environment. These models can be one-, two- or three-dimensional and can model systems at steady-state, where inputs and outputs are constant through time, or dynamically, where inputs and outputs are variable with time. The early measures of environmental system performance were deterministic predictions of ambient pollutant concentrations from environmental response models. However, these deterministic measures do not characterize all important aspects of system performance because they fail to account for uncertainty regarding the state of the prevailing environmental conditions and the uncertainty associated with the ability of the environmental response model to accurately model the system. Hashimoto et al. (1982) introduce alternative system performance indicators or indicators, such as reliability, vulnerability and resilience, that are risk-based and therefore can be used to account for the various inherent uncertainties associated with environmental response models. Reliability, vulnerability, and resilience characterize the frequency, magnitude and duration of the failure of water resource systems, respectively. Methods of determining the best or optimal management solutions are referred to as optimization or management models and include a variety of mathematical programming techniques and heuristic optimization methods. Environmental management models usually use environmental response models of the system being managed to predict the environmental response to the alternative management options. In many environmental systems, the existence of multiple management objectives requires the application of a multi-objective optimization technique to generate a variety of solutions that accurately and efficiently characterize the tradeoffs between the objectives. 1 Traditional solution techniques to water quality management problems require simplifying assumptions or rely on linear characteristics of the environmental response models. These traditional techniques and a variety of subsequent solution techniques are generally not applicable to case studies with combinations of (1) more complex nonlinear environmental system response models (e.g. dynamic models with two- or three-dimensions) (2) the quantitative treatment of uncertainty using risk-based system performance indicators or (3) multiple pollutants requiring simultaneous management. The growing popularity of artificial neural networks (ANNs) in water resources (see Maier and Dandy 2000) as environmental response models exemplify the fact that many water resource systems are nonlinear. ANNs are a nonlinear black-box modelling technique that Team' the relationships between input and output data and are developed after the functioning of neurons in the brain. Besides ANNs, many environmental response models are growing in size and scope and therefore demand more and more computational time to generate output predictions. Thomann (1998) believes that for effective ecosystem management the future of water quality response models lies in the incorporation of all component environmental response models into aggregate environmental or ecosystem response models that predict pollutant levels in all phases of the environment - water, air and land. In other words, the water quality response model would simply be a module that accepts inputs from, and provides outputs to, sediment, hydrodynamic, watershed and airshed transport and fate models. Almost certainly, such enormous aggregate response models will be highly nonlinear, three-dimensional, simulate numerous contaminants and be very computationally burdensome. The computational burden is exacerbated if these aggregate environmental response models are used within management models to evaluate the environmental response to many potential management options. Furthermore, when multiple management objectives exist, a number of candidate management alternatives that characterize the tradeoffs between the conflicting objectives are often generated. This is usually accomplished by repeatedly solving a number of management models over the range of management objectives. In this context, it is clear that parallel advancements are needed in the systems analysis area in order to use such complex, nonlinear and computationally expensive environmental response models within environmental management models. Robust, nonlinear, multi-objective optimization methods with the ability to incorporate system uncertainty characterizations and risk-based system performance indicators need to be further developed and then demonstrated as 2 applicable to complex case studies. Furthermore, due to the size of such environmental response models, the efficiency of any new optimization method will likely always have to be considered. This thesis demonstrates the use of Genetic Algorithms (GAs) in combination with reliability analysis techniques and a water quality response model for multi-objective water quality management under uncertainty. GAs are heuristic optimization techniques that are utilized here to search for optimal management solutions. The>reliability analysis techniques and water quality response model function together to estimate the risk-based system performance indicators by using either the First-Order Reliability Method (FORM) or Monte Carlo Simulation (MCS) for all candidate management solutions generated by the GA. The approach is demonstrated for carbonaceous biochemical oxygen demand (CBOD) - dissolved oxygen (DO) management on the Willamette River, Oregon, USA. The management objectives considered include the minimization of total CBOD waste treatment costs and the optimal attainment of the risk-based system performance indicators of reliability, vulnerability and resilience. The performance indicators are estimated with respect to meeting DO standards and the ambient DO is simulated using a QUAL2E water quality response model that has recently been developed for the Oregon Department of Environmental Quality (ODEQ) as part of the Willamette River Basin Water Quality Study (WRBWQS). The sources of uncertainty considered include the streamflow from four tributaries, river temperature and the water quality response model reaeration (Ka) and sediment oxygen demand (SOD) coefficients. This work, as well as the author's work in the complement papers (Maier et al. 1999; Tolson and Lence 2000; Vasquez et al. 2000) is significant and novel in a variety of ways: â€¢ FORM is used for the first time to estimate resilience and vulnerability â€¢ FORM and MCS estimates of reliability, vulnerability and resilience are thoroughly compared with respect to accuracy and efficiency. The methodology behind the efficiency comparison here is unique from other comparisons in the literature and the presentation of the advantages and disadvantages of FORM and MCS serve as general guidelines on the choice of the performance indicator estimation technique. â€¢ FORM and MCS performance indicator estimation techniques are combined with a GA to generate tradeoff curves for a variety of two-objective water quality management formulations. 3 â€¢ The Tradeoff Surface Representation (TSR) Algorithm is incorporated as an accurate and efficient multi-objective approach for approximating convex tradeoff surfaces. Comparison with the commonly used Constraint Method for multi-objective tradeoff curve development shows that the TSR Algorithm is superior since it produces a more representative tradeoff curve when the same number of tradeoff solutions are generated. The TSR Algorithm also provides a maximum amount of auxiliary information regarding the bounds of the tradeoff surface between identified tradeoff solutions. The scope of this thesis encompasses a number of important topics in water quality modelling and more generally, water resources. The goal of this work is to present viable methods and models for managing water quality systems under uncertainty. This work builds upon methods developed and applied in a number of fields including water quality management. Further, new techniques are developed and presented in this thesis that apply to all of these fields. The results of this work present a multi-objective framework for managing water resources, or even any system in general, under uncertainty. The balance of this thesis is structured as follows. Chapter 2 presents a broad review of the literature related to the various topics in this thesis. The methods of analysis, which are FORM, GAs, MCS and the TSR Algorithm, are described in detail in Chapter 3. Chapter 4 presents the Willamette River case study and outlines the data, details and assumptions related to the estimation of reliability, vulnerability and resilience. In addition, Chapter 4 also presents the methods used for the waste treatment cost data generation, the water quality management model formulations and the programming framework required to implement the management models. The results are presented and discussed in Chapter 5 in two sections. The first section compares the FORM-based estimates for the performance indicators with MCS estimates while the following section presents the water quality management model results and comparisons. Reliability and waste treatment costs are used as the objectives to compare the Constraint and TSR Methods and are also the objectives utilized in a FORM and MCS comparison for reliability estimation within the water quality management model. Chapter 6 summarizes the conclusions of this study and outlines future work that should continue in this area. 4 2 LITERATURE REVIEW The literature review presented here draws from research in a number of fields. In Section 2.1, a history of water quality management is presented. Next, Section 2.2 gives an introduction to the topic and sources of uncertainty in environmental response modelling. Section 2.3 is a review of the development and use of various risk-based system performance indicators in water resources. In Section 2.4, recent water resources applications of MCS and FORM for estimating the riskbased system performance indicators are described. GAs as an optimization tool in water resources is then reviewed in Section 2.5. In Section 2.6, some general and more recent advanced techniques in multi-objective optimization are discussed. Finally, in Section 2.7, the literature review concludes with a review of some of the more recent studies in water resources that present optimization models for managing water resource systems under uncertainty. 2.1 Water Quality Management Human induced contamination of natural water quality systems can have devastating environmental and even costly anthropogenic consequences such as algal blooms and fish kills. Point and non-point sources can be responsible for such events, but since the results of point sources are the most easily regulated sources, they tend to be the primary focus in water quality management. The general goal of water quality management is to design point source regulation programs for waterways or waterbodies such that a variety of societal objectives are maintained at optimal levels. Societal objectives can include total pollution abatement cost minimization, equity, administrative feasibility and maintenance of water quality at acceptable levels. The management decisions to be made are the level of waste treatment that each point source must provide before emitting the waste to the surface water body. The societal objectives are often combined together by treating one as the objective and the others as constraints. For example, the classical problem in surface water quality management is to make the set of decisions that will minimize total pollution abatement costs while still maintaining water quality to acceptable levels. Early solutions to this problem involved mathematical programming techniques where the cost is the objective and the limits on ambient water quality levels are constraints generated by a water 5 quality response model that describes ambient conditions as a function of the waste treatment levels at each point source. Initial work focused on deterministic problems and some of the earliest efforts used Linear Programming (LP) techniques linked with one-dimensional steadystate Streeter-Phelps DO response models to find least-cost solutions to wasteload allocation problems (Loucks et al. 1967; Revelle et al. 1968). The LP approach requires that all objectives and constraints be linear. Non-linear cost functions are usually approximated by piecewise linear curves while the linearity of water quality constraints results from the properties of the StreeterPhelps DO response model. The Streeter-Phelps model is based on superposition theory where all point source waste effluent effects are separate and additive. This characteristic allows the impacts of each point source on the ambient DO level to be described by linear transfer coefficients. Dynamic programming (DP) applications grew popular because they do not require that either the problem objective or constraints are linear (Liebman and Lynn 1966; Shih 1970). Subsequent work incorporates stochastic or uncertain elements of a water quality system using stochastic LPs (Lohani and Saleemi 1982) and stochastic DPs (Cardwell and Ellis 1993). In theory, optimal solutions to a variety of such problems can and have been found many times over. However, in practice these optimal solutions are almost always not implemented. Although there are many reasons for this, perhaps one of the most straightforward reasons is that in most countries, regulations and standards have been developed in ways that do not allow for optimal management strategies. Somlyody (1997) discusses the various impediments to implementing least-cost or optimal water quality management policies. Despite the lack of practical use, the continued research advancements and applications in this area are important because they repeatedly show economic incentives for rewriting regulations and they provide benchmark solutions with which all other management strategies can be compared. 2.2 Uncertainty in Environmental Response Modelling Environmental decision-making is quite often based on the results of complex environmental response models. Since these models are approximations or idealized representations of reality and most of their input and parameters are subject to various degrees of uncertainty, their predictions are also uncertain. Thus, effective environmental decision-making should be based on response model results that give an idea of the amount of uncertainty associated with the response model's predictions rather than a simple deterministic result. Indeed, much of today's 6 modelling efforts attempt to quantify the amount of uncertainty associated with environmental response model predictions. Before such model prediction uncertainties can be quantified however, the sources of uncertainty must be determined and characterized since the method of quantification depends on the source of uncertainty (Morgan and Henrion 1990). A number of studies define various sources of uncertainty that arise in environmental modelling and management decisions. Burges and Lettenmaier (1975) divide modelling error into two types of errors in a water resources context. Type I error results from the use of an inappropriate model within proper deterministic parameters and can be further subdivided into modelling error associated with the incorporation of many processes into one and the error associated with inappropriate model selection. Type II error is a result of the use of uncertain parameters in a model that perfectly characterizes the system. Tung and Hawthorn (1988) give a similar characterization of uncertainty in the context of water quality modelling where the sources include inherent, model and parameter uncertainty. Inherent uncertainty refers to the random attributes of the physical, biological and chemical properties of the pollutant transport process through space and time. Model uncertainty is a result of using an idealized model to represent a more complex process. Total parameter uncertainty results from parameter estimations based on limited field data and even imprecision and bias in the measurements of parameters. A thorough review by Beck (1987) also presents problem areas in uncertainty analysis within water quality modelling and summarizes other attempts to classify uncertainty sources. He concludes that larger or more comprehensive water quality models are likely to produce highly uncertain predictions of future ambient conditions. A complete treatment of quantitatively accounting for uncertainty in decision-making is given in Morgan and Henrion (1990). They first categorize the types of modelling quantities that are defined in policy models (which include environmental response and management models) and then advocate that the only modelling quantities whose uncertainty should be expressed by probability distributions is that of empirical quantities. Empirical quantities are model parameter or input values that, in principle, can be measured at some point in time and examples include reaction rates and natural inputs such as temperature and flow. Morgan and Henrion (1990) further classify the sources of uncertainty in empirical quantities into categories such as statistical variation (random and systematic measurement errors), disagreement and variability. Disagreement refers to the fact that different estimates or measures of the same quantity can 7 result when different experts provide the estimate or measure. Variability is further subdivided into the variability of quantities over time and space and the empirical uncertainty in quantities due to incomplete scientific data. This distinction is further explained by considering that the variability of a quantity can be described by a probability distribution where as an example of the empirical uncertainty affiliated with that quantity would be our uncertainty about the specific probability distribution describing the variability of that quantity. Furthermore, while variability is a source of uncertainty in empirical quantities, it is fundamentally different than empirical uncertainty and thus the two should be distinguished and accounted for differently. Morgan and Henrion (1990) also discuss the uncertainty associated with model form. They state that, except for special cases, uncertainty about model form is a result of conflicting beliefs among experts about the underlying structures and assumptions in the model. Therefore, they advocate that such uncertainty is better investigated parametrically rather than probabilistically via probability distributions. The characterization of empirical quantities with probability distributions can be utilized to generate a probability distribution describing the model output. This probabilistic output can be used to describe the system performance under uncertainty. 2.3 Risk-Based System Performance Indicators in Water Resources The first attempts to describe the performance of water resource systems under uncertainty focused on characterizing the mean and variance of the system output predictions. Hashimoto et al. (1982) show that there are a number of other important aspects of system performance under uncertainty. Since most water resource and environmental systems can be described as being in either a satisfactory or unsatisfactory (e.g. system failure) state, system performance measured in terms of failure or non-failure are more relevant to decision-makers. Hashimoto et al. (1982) first argue that attention should be focused more on the system performance during failure rather than on the system performance during non-failure states. They provide a further simple example to show that two systems can have the exact same mean and variance values of system performance (e.g. ambient DO concentrations) over the same time period but also have a different number of system failures (e.g. failure to meet the DO standard). 8 Hashimoto et al. (1982) introduce and operationalize the risk-based system performance indicators of reliability, vulnerability and resilience and illustrate their use with a water supply reservoir under various operating policies. Reliability is the probability that a system is in a satisfactory or non-failure state and is the complement to the probability of being in a failure state. Vulnerability is the likely magnitude or significance of failure if one occurs. Resilience is equivalently defined as (i) the inverse of the expected length of time a system's output remains in a failure state or (ii) the conditional probability that the system will recover from failure in a single time step, given that failure occurred in the previous time step. Similar concepts have been previously introduced by Fiering (1969; 1976) to assess water supply systems and assist in water resources planning and by Holling (1973) to describe ecological resilience. Holling (1973) describes ecologically resilient systems as having the ability to avoid being shifted to a less preferred or non-equilibrium state when experiencing ecological stress or variability. Duckstein and Bernier (1986) present a detailed list of seven engineering performance indices related to risk in water resources. These include the above risk-based performance indicators, which are defined similarly to those by Hashimoto et al.(1982), as well as other economic indices. The risk-based system performance indicators given by Hashimoto et al. (1982), or variations thereof (Moy et al. 1986; Burn et al. 1991), have been used to assess reservoir operating policies (Hashimoto et al. 1982; Moy et al. 1986; Burn et al. 1991; Kumar et al. 1996), to assess the performance of water distribution systems (Zongxue et al. 1998) and to characterize regional droughts (Correia et al. 1986). Applications in water quality for reliability estimates are more abundant than research quantifying vulnerability or resilience. Any study that characterizes the risk of failing to meet a specified water quality standard is equivalent to one that characterizes the complement reliability. Recent studies using more than just reliability as a system performance measure include the author's work in Maier et al. (1999) where the risk-based performance indicators are estimated over varying CBOD wasteloads with respect to meeting a DO standard (this work is also presented here in Chapter 5) and the work of Chan-Yan (2000) who estimates the reliability and resilience of meeting a turbidity standard during severe drawdown conditions in a water supply reservoir. Logan (1990) evaluated the reliability, vulnerability and resilience, as defined by Hashimoto et al. (1982), of several rivers in Ontario, Canada with respect to contaminants such as total phosphorous, nitrates, copper and suspended solids. The 20-year monthly history of record is 9 used to derive reliability, vulnerability and resilience curves over various threshold contaminant levels. These curves show that the performance indicators are not conflicting over various threshold levels (i.e. improving one measure improves the other measures and all of the best individual performance measure values occur at the same threshold level). Logan (1990) also modified the resilience definition by introducing what is essentially a conditional resilience estimate where the resilience estimate is conditional on the length of the failure period rather than being conditional on the system being in failure only during the previous time step. He demonstrates this by showing that for two different threshold levels the resilience can either increase or remain essentially constant as the length of failure increases. Most of the studies discussed above use simulation to estimate the various performance indicators (with the exception of Logan (1990)). The data sets used in the simulations are either relatively short historical records or synthetically generated data. A more recent approach developed and applied in the author's work in Maier et al. (1999) and also applied by Chan-Yan (2000) is to use reliability analysis methodology and the First-Order Reliability Method to estimate reliability, vulnerability and resilience. 2.4 First-Order Reliability Method (FORM) in Water Resources FORM is a reliability analysis technique used to assess the probability that system performance, defined in terms of its load and resistance, is in a non-failure state. FORM was originally developed to assess the reliability of structures (Hasofer and Lind 1974; Rackwitz 1976). However, use of the load-resistance analogy for water resources problems has been discussed by a number of authors, including Duckstein and Bernier (1986) and Kundzewicz (1989). For example, in the water supply case, water demand corresponds to system load and supply capacity to system resistance, whereas in the water quality case, pollution load and a particular water quality standard correspond to the system's load and resistance, respectively. The FORM technique is outlined in Section 3.1 while other detailed descriptions can be found in Madsen et al. (1986), Sitar et al. (Sitar et al. 1987), Melching (1992) and Skaggs and Barry (1997). Other reliability analysis techniques include MCS, Mean-value First-Order Second-Moment analysis (MFOSM) and the Second-Order Reliability Method (SORM). Detailed descriptions of the MFOSM approach, as well as comparative studies between MFOSM and FORM, are given 10 by Tung (1990) and Melching and Anmangandla (1992) while a detailed description of SORM is given by Madsen et al. (1986). FORM applications in water resources engineering are primarily to groundwater problems for analyzing subsurface flow and contaminant transport (Sitar et al. 1987; Jang et al. 1994; Skaggs and Barry 1997), although there have been some surface water applications. For example, Tung (1990) compares the performance of FORM and MCS for evaluating the probability of violating various DO standards for a hypothetical case study. A similar study, carried out by Melching and Anmangandla (1992), uses the hypothetical DO case studies of Burges and Lettenmaier (1975) and Tung and Hathhorn (1988). In both papers, the performance of FORM is very similar to that of MCS. Melching et al. (1990) uses FORM to determine the uncertainty of the peak discharge predictions obtained from a rainfall-runoff model for the Vermillion River watershed, Illinois, USA. Melching (1992) carries out a comparison between FORM and MCS for the same case study and finds good agreement between FORM and MCS predicted failure probabilities for a wide range of storm magnitudes and types. These surface water applications of FORM only estimate the system reliability and assess the FORM predictions against MCS. Recent work by the author in Maier et al. (1999) pioneers the use of FORM for the estimation of vulnerability and resilience. To the author's best knowledge, there are only three studies that utilize FORM for reliability estimation within an optimization model in water resources. The author's work in Vasquez et al. (2000) and Tolson and Lence (2000) uses FORM as the reliability analysis technique within multi-objective water quality management models to find the cost-reliability tradeoff and the CBOD total wasteload-reliability tradeoff respectively, for case studies on the Willamette River in Oregon, U.S.A. Xu and Goulter (1999) use FORM to estimate reliability in a least-cost design of water distribution networks. 2.5 Genetic Algorithms (GAs) in Water Resources GAs are heuristic optimization algorithms that are based on Darwinian evolution and survival of the fittest in the natural environment. GAs were first developed and introduced by Holland (1975) while a more recent popular introductory text on GAs is Goldberg (1989). They are an unconstrained optimization technique that requires constraints to be accounted for by functions that penalize solutions that violate the constraints. These functions are called penalty functions. 11 GAs have been applied successfully to optimization problems across a large number of disciplines including water resources. They appear to be one of the most common optimization methods in water resources given their preponderance in the recent body of literature. Generally speaking, GAs have been applied to almost every type of optimization problem encountered in water resources. For example, they have been used successfully as an optimization technique in water quality, rainfall-runoff and groundwater model calibration (Wang 1997; Mulligan and Brown 1998; Solomatine et al. 1999), pipe network design and rehabilitation (Goldberg and Kuo 1987; Halhal et al. 1997), irrigation system management (Chen 1997) and reservoir management (Oliveira and Loucks 1997; Wardlaw and Sharif 1999). Example applications of GAs for groundwater management are among the most common applications (see Ritzel et al. 1994; Cieniawski et al. 1995; Huang and Mayer 1997; Wang and Zheng 1997; Aly and Peralta 1999a; Aly and Peralta 1999b; Morshed and Kaluarachchi 2000). These applications include the optimization of groundwater well pumping rates, and sometimes the well locations, for single and multiple objective problems in either deterministic or uncertain systems. The work by Aly and Peralta (1999b) analyzes optimal aquifer cleanup strategies under uncertain parameter values and will be discussed in more detail in Section 2.7. GA applications in surface water quality management are less numerous. Burn and Yulianti (1999) demonstrate the application of GAs to a surface water quality management problem and evaluate the tradeoff between total treatment cost, water quality levels and equity for a set of steady-state critical stream water quality conditions. Chen and Chang (1998) use a GA in combination with fuzzy set theory in a multi-objective framework to minimize biochemical oxygen demand (BOD) waste treatment costs, maximize the assimilative capacity of the river and maximize the economic value of the river flow. GAs are used as the optimization technique in the Vasquez et al. (2000) and Tolson and Lence (2000) water quality management studies mentioned in Section 2.4. 2.6 Multi-Objective Optimization Multi-objective problems arise where systems are managed for two or more conflicting objectives. Two objectives are conflicting if the continuous improvement of one objective eventually results in the degradation of the other objective. The basic idea in managing a multi12 objective system is to find a solution, or set of decision variables, that gives the 'best' of all the objectives. The choice of the best solution is usually determined by the decision-maker, whereas the generation of alternatives, from which to select the best solution, is usually the task of the systems analyst. At this point it is necessary to introduce some simple multi-objective terminology. A solution in which no single objective can be improved without degrading any other objective is termed a Pareto-optimal solution. The Pareto-optimal solutions are the set of solutions from which the decision-maker selects his/her preferred management strategy. A non-Pareto-optimal solution would not be selected as the preferred management strategy since at least one solution would exist in which the value of one objective could be improved without degrading the value of any of the other objectives. The Pareto-optimal surface, or tradeoff surface, is the surface of connected Pareto-optimal solutions in objective space. An ideal point in objective space is one in which all objectives are optimized individually and is ideal because, while it is the best individual set of objective values, no unique solution will exist such that the ideal point is simultaneously reached for all objectives. This section reviews the common methods of multi-objective optimization and more recent novel multi-objective approaches. Traditional approaches to multi-objective optimization all require that multiple objectives be combined into a single objective function in order to generate a single Pareto-optimal solution. In order to generate a set of Pareto-optimal solutions and construct an approximate tradeoff surface, this process must be repeated for a series of different single objective functions. Narayanan and Azarm (1999) identify three general approaches to multiobjective problems. Posteriori methods involve the generation of the set (or representative subset) of Pareto-optimal solutions for presentation to the decision-maker who then selects his/her preferred solution. Interactive methods imply that the decision-maker and systems analyst work together to identify the preferred solution without identifying a completely representative set of Pareto-optimal solutions. Priori methods differ from posteriori and interactive methods in that instead of generating multiple Pareto-optimal solutions, they require the decision-maker to express his/her preferences in terms of a single combined objective function prior to optimization. Since this research focuses on the identification of multiple Pareto-optimal solutions, only posteriori and a few interactive methods will be discussed further. 13 Eschenauer et al. (1990) provides a review of many of the basic multi-objective concepts and a variety of more advanced interactive multi-objective approaches. One common method to generate a range of points on the tradeoff surface is the Weighting Method where the objectives are combined to a single objective function that is a weighted sum of the objectives. This is repeated with a number of weighting sets to generate a range of points on the tradeoff surface. However, this method only produces a range of Pareto-optimal points when the tradeoff surface is convex and Das and Dennis (1997) show that even for convex problems a smooth spread of weights frequently fails to produce an even spread of points on the tradeoff surface. Another common technique that overcomes these weaknesses is the Constraint Method. The Constraint Method requires constraining m-l of the m objectives to meet a given threshold and then optimizing the unconstrained objective. This process is repeated at equal intervals for various thresholds of the constrained objectives to generate a set of points on the tradeoff surface. A more detailed explanation of the Weighting and Constraint Methods can be found in Cohon(1978). Two studies proposing interactive approaches are also noteworthy. Palli et al. (1998) develop the Multistage s-Inequality Constraint Method for effective interactive multi-objective decisionmaking. The difference between this and the Constraint Method is that this method is applied in stages where, at the end of each stage, the decision-maker essentially selects the solution that lies in the area of most interest and then the process is repeated with a smaller range of constraint values centered around the selected solution, so as to iteratively locate the preferred solution area. Another difference is that the constraints are spaced randomly rather than evenly for each constrained objective. Fandel (1972), as described Eschenauer et al. (1990), (the original work by Fandel is contained in a manuscript written in German) introduces a unique interactive approach to multi-objective decision-making. One iteration of this approach involves first finding the tradeoff surface point between m other tradeoff points (Pi, P2, to the shifted hyperplane defined by the points (Pi, P , 2 P,â€ž) that is tangent P,â€ž). Then, the decision-maker defines a reduced set of bounds on the objectives that contains the tangent point on the tradeoff surface. The intersection points of these bounds with the tradeoff surface are then determined and used to define another hyperplane and then find the tradeoff point that is tangent to it. This process is repeated until the decision-maker selects a preferred solution. The determination of the intersection points and the tangent points are not discussed beyond their graphical equivalent in 14 Eschenauer et al. (1990). However, the determination of each point would require the solution of separate single objective functions. Das and Dennis (1998) present a detailed account of an innovative nonlinear posteriori method to generate an even spread of points along convex tradeoff surfaces called the Normal-Boundary Intersection (NBI) method. The NBI Method works by first defining a hyperplane with the bounds of the objectives and then defining regularly spaced normal vectors along the hyperplane in the direction of the ideal point. Then, for each vector, a nonlinear program (NLP) finds the solution where the vector and tradeoff surface intersect. This is achieved by constraining the NLP to maximize the distance along the normal in the direction of the ideal point. The method is shown to produce an even spread of points along the tradeoff surface independent of objective function scales, works equally well for convex and non-convex curves and can be easily extended to any number of objectives. Inherent in the NBI formulation according to Das (1999) is the ability to find the 'knee' or the 'bulge' in an w-objective tradeoff surface. The knee is the Pareto optimal point that is furthest from the hyperplane defined by the bounds of the objectives in the direction of the ideal point. Based on the NBI formulation, the knee is found by allowing the point of origin of a normal vector to be an optimization variable instead of being fixed. Das (1999) shows that this method can be used for convex tradeoff surfaces. Besides the above multi-objective methods that reduce all objectives to a series of single objective functions there have been recent enhancements to GAs for multi-objective optimization that generate the entire surface without constructing a series of single objective functions. These GAs are referred to as Multi-Objective Genetic Algorithms (MOGAs). A variety of types exist but the basic idea is that the MOGA attempts to identify the entire tradeoff surface in one optimization run instead of generating single points on the tradeoff surface. Goldberg (1989) first proposed the Pareto-based approach for MOGAs where the fitness or performance of a solution is dependent on whether that solution is Pareto optimal. Zitzler and Thiele (1999) describe and compare the performance of five variations of multi-objective evolutionary algorithms. 15 2.7 Recent Management Models for Water Resource Systems Under Uncertainty In the study by Chen and Chang (1998) previously described in Section 2.5, a GA is used in a multi-objective water quality management model. The tributaries are treated as large point sources and the levels of waste treatment for each tributary are the decision variables. They use a weighting technique to combine all objectives for the multi-objective approach and fuzzy set theory to consider the uncertainty due to the vagueness related to describing the water quality and pollutant abatement goals. This study however, does not consider any other sources of uncertainty and treats the uncertainty as fuzziness rather than probabilistically. Takyi and Lence (1999) summarize three of the most commonly used methods for incorporating uncertainty into environmental management models and thus generating tradeoff curves for the objectives of total waste treatment cost and the reliability of meeting water quality goals. These methods include Chance-Constrained Programming (CCP) (Lohani and Thanh 1978), Combined Simulation-Optimization (CSO) (Fuessle 1987) and multiple realization- (MR-) based approaches (Wagner and Gorelick 1989; Morgan et al. 1993). CCP accounts for random variables in the problem by incorporating them in probabilistic constraints and then solving the resulting mathematical program. The CSO approach consists of generating a number of possible scenarios and then for each scenario, optimizing the resulting deterministic system. The MR- based approach treats each of a number of realizations for the random variables as a deterministic constraint in which the solution must satisfy the constraint (i.e. a water quality standard). To generate a number of Pareto-optimal solutions at different reliability levels, critical realizations are identified and then removed and the problem is resolved with a reduced set of constraints in an iterative manner. While the CCP and MR approaches can explicitly determine the treatment cost-reliability tradeoff curve, neither approach is efficient for solving nonlinear problems, and CCP has the additional disadvantage of requiring simplifying assumptions about the input information and its distributions. The CSO approach accommodates nonlinear problems, but does not explicitly generate the treatment cost-reliability tradeoff, and may not produce cost efficient decisions at certain reliability levels (Fuessle 1987). Another approach to optimizing uncertain systems is to embed a MCS of the system (or any other reliability estimation technique such as FORM) in an optimization framework. In other words, this approach explicitly estimates the reliability and cost for each solution. This approach 16 is most closely related with the MR-based approach. In order to compare the approaches a few relevant studies are examined in more detail. One of the first studies to employ MRs is by Wagner and Gorelick (1989) who used it to design reliable aquifer remediation design under spatially variable hydraulic conductivity. Their focus is on reliable design only and as such they develop single reliable optimal solutions instead of a treatment cost-reliability tradeoff. They conclude that using a set of only 30 realizations of the uncertain parameters and constraining the solution to satisfy the cleanup requirements for all 30 results in very reliable designs. This conclusion is based on checks of the optimal solutions' reliability via a post-optimization MCS with 1000 realizations. Morgan et al. (1993) develop an extended approach to that of Wagner and Gorelick (1989) for a similar type of aquifer remediation problem under hydraulic parameter uncertainty, where the extended approach functions to develop tradeoff curve points over a range of reliabilities. Multiple tradeoff points â€¢are generated by successively solving the optimization problem with a reduced set of constraints (e.g. a reduced number of realizations) and then assigning a nominal reliability to these successive solutions. The most critical realization must be identified as the one to be removed during the next stage and the nominal reliability is then the proportion of the number of present model constraints over the original number of realizations used. Takyi and Lence (1999) improve on this method in a least-cost CBOD wasteload water quality framework by incorporating new methods to select the various realizations to be removed at each stage, such as using an ANN for identification of important or critical realizations. The increased complexity in the method used by Morgan et al. (1993) and Takyi and Lence (1999) caused by having to identify critical realizations to remove at each stage is somewhat undesirable. Further, the nominal reliability estimates for the successive optimal solutions are not necessarily the same as they would be if all original realizations were used to evaluate the reliability of the optimal solutions. In other words, the nominal reliability is in fact an estimate of the MCS estimate of reliability using the entire set of realizations. Lastly, with these methods it is unclear how to analyze problems in which there is no feasible solution that will achieve a reliability of 1.0. For example, it is easy to imagine that in many environmental systems, due to background conditions and natural variation, even the most expensive solution (i.e. our best effort at cleaning up the system) will not result in a solution that will satisfy all constraints on the realizations. It is also not clear whether MR-based approaches can be extended to more than two 17 objectives to include objectives such as vulnerability or resilience. While the general structure of the method makes it possible to include vulnerability or resilience in place of reliability, their inclusion may complicate the identification of critical realizations. A recent paper by Aly and Peralta (1999b) on the optimal design of aquifer cleanup systems attempts to improve on the basic MR approach used by Wagner and Gorelick (1989). The improvement is achieved by using a GA as the optimization technique and approximating the groundwater flow and contaminant transport model with an ANN to approximate the complex system response model by way of a response surface. The decision variables for this problem are the discrete pumping rates and possible well locations and the objectives are to maximize the cleanup strategy reliability and minimize the total pumping rates. Their approach reduces the number of realizations used and is reformulated to minimize the worst level of contamination observed from the set of realizations subject to various constraints on the total pumping rates. This reformulation is analogous to maximizing the reliability and using the Constraint Method to generate a representation of the reliability-total pumping rate tradeoff curve. Although this approach does embed a small MCS within the optimization model, the reliability of each optimal solution is not explicitly estimated within the optimization model. The ANN model predicts the worst level of contamination for each realization and varying the constraint level for the total pumping rate produces new points on the reliability-total pumping rate tradeoff curve. Actual reliability however, is estimated post-optimization for each optimal solution using 250 to 300 MCS realizations. For this problem, a small number of realizations (10) relative to other studies are adequate to determine optimal pumping strategies. The study by Aly and Peralta (1999b) does not generate the reliability-total pumping rate tradeoff curve in the same manner as Morgan et al. (1993) and Takyi and Lence (1999). The selected constraint level on the total pumping rate determines the solution reliability level. Therefore, a trial and error method must be used to determine the appropriate constraint levels that result in a range of reliability values that are of interest. In this study, this does not pose a significant problem. However, in studies where the relevant constraint levels and the corresponding reliabilities are less well known, it may become more difficult and time consuming to find the relevant constraint levels by trial and error. 18 Work by the author in Vasquez et al. (2000) avoids some of the problems encountered in the MR approaches above because GAs are used in place of mathematical programming techniques and a large number of MCS realizations or FORM estimates the reliability of each solution within the optimization model. Vasquez et al. (2000) combine a GA with FORM and MCS to develop a total CBOD waste treatment cost-reliability tradeoff curve. They approach this multi-objective problem with the Constraint Method, where reliabilities are constrained at minimum levels at regular intervals and the total waste treatment costs of ten point sources are minimized, to generate this tradeoff. They also show that although the GA-FORM tradeoff curve is very similar to that of the GA-MCS tradeoff when 5000 MCS realizations are used for each reliability estimate, the GA-FORM method is much more efficient. The use of a simple Streeter-Phelps DO model made 5000 realizations computationally tractable for this problem. Tolson and Lence (2000) use a GA-FORM based model with a more complex water quality model than the one used in Vasquez et al. (2000) to develop the total CBOD wasteload-reliability tradeoff curve. Maximizing the total wasteload is used as a surrogate for minimizing waste treatment costs. It is shown that similar tradeoff curves result when total wasteloads and reliabilities are interchanged as the constraint and objective. However, the formulation of maximizing reliability with total wasteload as constraints is found to be more efficient because this formulation allows the model to attempt to satisfy the wasteload constraints implicitly rather than relying heavily on a penalty function. 19 3 METHODS OF ANALYSIS 3.1 Reliability Analysis The performance of any engineered system can be expressed in terms of its load (demand) and resistance (capacity). If X = {X\, Xi, X) T n is the vector of random variable values that influences a system's load (L) and/or resistance (R), the performance function, G(X), is commonly written as: G(X) = R-L (1) In some water quality cases such as DO, a low load value leads to system failure. Thus, in the most general case, the performance function must be defined such that values below zero characterize system failure. The failure (limit state) surface, G = 0, separates all combinations of X that lie in the failure domain (F) from those in the survival domain (S). Consequently, the probability of failure, pf, is given as (Sitar et al. 1987): (2) Pf G(X) <0 where fx(x) is the joint probability density function (PDF) of X. In most realistic applications, the integral in Equation (2) is difficult to compute (Sitar et al. 1987). Approximate solutions can be obtained by using a variety of techniques including MCS, MFOSM, FORM, and SORM. The most accurate reliability estimation method is MCS with a very large number of realizations. This research focuses on the comparison of FORM and MCS for the prediction of reliability, vulnerability and resilience. 20 3.1.1 FORM As mentioned in Section 3.1, the objective of FORM is to obtain an estimate of the integral in Equation (2), and hence the probability of failure. The basic idea is to compute the "reliability index", 6, which is then used to obtain the probability of failure using: Pf=<K-P) (3) where cP( ) is the standard normal cumulative distribution function (CDF). In the rc-dimensional space of n random variables, B can be interpreted as the minimum distance between the point defined by the values of the n variable .means (mean point) and the failure surface (Figure 3.1). Consequently, B may be thought of as a safety margin, as it indicates how far the system is from failure when it is in its mean state. The point on the failure surface closest to the mean point is generally referred to as the design point, X (e.g. Figure 3.1), which may be thought of as the most likely failure point. In other words, the design point yields the highest risk of failure among all points on the failure surface. Determination of the design point, and hence 8, is a constrained nonlinear minimization problem (Tung 1990). Suitable optimization techniques include the Rackwitz-Fiessler Method (Madsen et al. 1986), the Generalized Reduced Gradient Algorithm (Cheng 1982) and the Lagrange Multiplier Method (Shinozuka 1983). In this thesis, FORM is implemented using a slightly modified version of a general RELiability ANalysis software package called RELAN (Foschi and Folz 1990), which uses the Rackwitz-Fiessler Algorithm. Equation (3) is only exact if (i) the elements of X are uncorrected normal variables with a mean of zero and a standard deviation of one and (ii) the failure surface is a hyperplane. These conditions are rarely met in realistic applications. The approach taken in FORM to deal with the first of these problems is to transform all random variables (X\, X , 2 uncorrelated standard normal variables (V\, Vi, X ) to the space of n Vâ€ž). In the original RELAN version, the method of Der Kiureghian and Liu (1986) is used to perform this transformation, as it accounts for the correlation structure between the random variables. However, in this study, a modified version of RELAN is used where the random variable correlations remain unchanged during this 21 transformation. The reasons for this will be discussed further in Section 4.2.1.4. The second condition cannot be accounted for exactly by FORM and the failure surface is approximated by its tangent hyperplane at the design point in standard normal space, V*, using First-Order Taylor Series Expansion (Figure 3.1). Consequently, the probability of failure obtained using FORM is only an approximation, unless the performance function is linear. The degree of nonlinearity in the performance function, and hence the accuracy of FORM, is problem dependent. first order approximation of failure surface design point Figure 3.1. F O R M approximation of the failure surface in standard normal space The sensitivity of the probability of failure estimate obtained with FORM to the uncertainty in the random variables can be determined easily by computing the direction cosines of the random variables in the standard normal space. For example, the direction cosine of the first random variable is defined as: cos # i = V *I j3 x (4) where 0\ is the angle between the vector of the design point, V, and the axis of the transformed random variable V\ and V\ * is the value of V\ at the design point. The sum of the squares of the direction cosines for all of the random variables is equal to one. The square of the direction cosine for a given random variable provides information about the sensitivity of the reliability to changes in that variable near the design point. The significance of the given random variable increases with increasing values of the square of the direction cosine. 22 3.2 FORM and Monte Carlo Simulation (MCS) Estimators of Reliability, Vulnerability and Resilience This section presents the operational definitions of the performance indicators, reliability, vulnerability and resilience, and discusses means of evaluating them using FORM and MCS. The FORM operational definitions are developed in the author's work in Maier et al. (1999). Also included is a discussion on the computational efficiency and accuracy of FORM and MCS for evaluating each measure. 3.2.1 Reliability Reliability is a measure of the probability of system survival. Hashimoto et al. (1982) define the reliability of a system, Â«, at time t as: a = ?r{X, eS} (5) which is the complement of the probability of failure. Therefore, the term 'probability of failure' is sometimes used throughout this thesis where conversely, reliability could also be used. Using Equation (3), reliability can thus be estimated as: a=\-p =l-C\-/3) f = aXJ3) (6) MCS can also be used to approximate the integral in Equation (2) by repeatedly sampling values of the random variables and evaluating the performance function in Equation (1). The reliability by MCS is then the ratio of the number of samples where G(X) > 0 over the total number of MCS samples taken. 3.2.2 Vulnerability Vulnerability is a measure of the magnitude of a system's failure. Hashimoto et al. (1982) define vulnerability, v, as follows: 23 (7) jeD where e,- is the probability that the system performance variable, L, is in discrete failure state j and Wj is a numerical indicator of the severity of failure state j (see Figure 3.1). If the discrete failure states are bounded by a hierarchy of failure levels, i?i < R < R3 < ... < Rj, e is given by: 2 ej = ?r{Rj <L<R } 7 = D (R ) -D (Rj) J+l L J+1 (8) L where D is the CDF of the load, L. Consequently, vulnerability is given by: L v = .(D (R^)-D (Rj)) L (9) L Using FORM, the probability of failure is given by: p = QX-B) = Pr{G(X ) < 0} = Pr{Z > R} = 1 - Pr{Z < R} = 1 - D (R) f t L (10) Consequently: D (R)=\-^rP) L = ^j3) (11) Combining equations (9) and (10), vulnerability can be expressed in terms of the reliability index, B, as follows: v = Â£ W j (0(/? ) -<P(fl)) = W j ) -^(-^ )) J+1 J+l jeF (12) jeF where ^ is the reliability index for resistance level Rj. Figure 3.1 presents the vulnerability definition in Equation (12) in a graphical manner. As pointed out by Melching et al. (1990) and Skaggs and Barry (1997), Equation (11) can also be used to obtain points on the CDF of the system performance variable, L, by repeating the reliability analysis for a range of values of system resistance, R. 24 MCS can also be used to estimate vulnerability by recording the system load that results from each MCS realization of the random variables and then counting the number of times the system fails with respect to the various resistance levels of interest. Each count divided by the total number of MCS samples gives the probability of failure with respect to each resistance level. Vulnerability is then calculated directly with Equation (9). 0 /?, R 2 R, Load(L) Figure 3.2. Cumulative distribution function (CDF) of load (Â£) and schematic representation of multiple failure states for vulnerability estimation 3.2.3 Resilience In water resources engineering, resilience has generally been used as a measure of how quickly a system recovers from failure, once failure has occurred. Hashimoto et al. (1982) give two equivalent definitions of resilience, y. One is a function of the expected value, (E[ ]), of the length of time a system's output remains unsatisfactory after a failure, 7}, as in Equation (13). The other is based on the probability that the system will recover from failure in a single time step (Equation (14)). 25 r = 1 (13) E[T ] f y = ?v{X M e S X t e F) = ?r{X e F g/K? e S} Pr{X e F} t (14) t In water resources applications, Equation (13) has generally been used to obtain estimates of resilience. This is undertaken by examining a time series (real or synthetic) of the system performance variable and counting the number of consecutive time steps the system remains in failure, once failure has occurred. FORM can be used to obtain estimates of resilience based on the conditional probability definition of resilience in Equation (14), as outlined below. In many instances in structural engineering, there is a need to consider multiple failure modes. For example, a beam may fail in bending or in shear and a retaining wall may fail by overturning or by sliding. If there are two failure modes, the total probability of failure is given by: P/ = PP+Pfl-Pfii (15) = Pr{G, < 0} + Pr{G < 0} - Pr{Gi < 0 and G < 0} 2 2 where pj\ and pp_ are the probabilities of failure due to failure modes 1 and 2, respectively; pf\ is 2 the joint probability of failure for failure modes 1 and 2 and G\ = G(X\) and G = G(X ) are the 2 2 performance functions for failure modes 1 and 2, respectively. The failure probabilities for the individual failure modes (p/i and pp) can be obtained using Equation (3). The joint probability of failure, pj\ , is given by Madsen et al. (1986) as: 2 Pa pn f $(-/?â€ž-/? ;A ) = 2 2 = $(-A)0(-A) + (16) o where @( , ;p) is the CDF for a bivariate normal vector with zero mean values, unit variances and correlation coefficient p and (//( , \p) is the corresponding PDF. The integral in Equation (16) is generally obtained numerically. The correlation coefficient needed to evaluate this integral, p\ , is calculated using Madsen et al. (1986): 2 26 where V\ and V are the design points in standard normal space for failure modes 1 and 2, 2 respectively. If the performance functions are defined as: G\=R,-L, (1 ) 8 (19) G = Li+i - R +i 2 t the corresponding individual and joint failure probabilities respectively, are given by: Pfl=Vr{XteF} = $(-B ) (20) l p = Pr{X, i e S} = <8(-/fc) (21) + fl p 2 = Pr{X, e FandX fl t+l e S} = Q\.p u -fc ) (22) pn It should be noted that the conventional definition of the performance function (see Equation (1)) is used in Equation (18). However, the order of L and R is reversed in Equation (19), so that the probability of failure, as defined in Equation (2), is actually the probability that the system will return to a non-failure state (see Equation (21)). Combining Equations (14), (20) and (22), resilience is given by: = $(-#,-A; A ) 2 ( 2 3 ) if system load and resistance are stationary processes, L, â€” L \, R = R i and pj\ = 1 - pp. t+ t l+ Persistence in the time series is accounted for by the lag one auto- and cross-correlations 27 between the elements of X and X . t t+l As presented in Hashimoto et al. (1982), if X and X +i are t t statistically independent, resilience is equivalent to reliability, and is given by: <K-A) r <*-A) ~ < 2 4 ) MCS can also be easily used to evaluate system resilience in accordance with Equation (14). This requires that each MCS resilience realization be composed of two consecutive sets of random variables corresponding to two consecutive time steps. The generation of random variables in two consecutive time steps is considered a single event. The random variables associated with the first time step are first used to evaluate the performance function for time step one in Equation (18) and the result is recorded. If the system fails in time step one, the variables associated with time step two are then used to evaluate a second performance function with the same form as Equation (18) and the result is also recorded. If the system does not fail in time step one, the performance function in the second time step is not evaluated and another MCS resilience realization is evaluated. The number of system failures in time step one divided by the total number of MCS resilience realizations can give the probability of failure in time step one (the denominator in Equation (14)). However, the resilience estimate, or the conditional probability of success in time step two given failure in the first time step, can be calculated directly by the number of successes in time step two divided by the number of failures in time step one. 3.2.4 Accuracy and Efficiency of FORM and MCS Estimates of the Performance Indicators As mentioned in Section 3.1.1, the FORM-based reliability estimate is rarely exact and the accuracy decreases as the nonlinearity of the failure surface increases. Thus, in all applications, FORM reliability estimates are based on an approximation where the accuracy of the approximation is problem dependent and can only be assessed by comparisons with reliability estimates obtained by MCS. In contrast, the efficiency of the FORM estimate of reliability can be determined by observing the convergence rate of FORM. The approximate number of performance function evaluations required in FORM as implemented by RELAN for each reliability estimate (7vÂ» is given by the following relationship: 28 N = (2Â«+l) + q(2n+\) = (q+l)(2n+l) (25) F where the first (2Â«+l) term is the minimum number of performance function evaluations required by an initialization procedure used to find good starting points for FORM; n is the number of random variables in the problem and q is the number of FORM iterations required to find the minimum B. Each FORM iteration requires 2n+l performance function evaluations because the performance function is evaluated for initial calculation and the gradient is calculated twice for each random variable. The minimum number of FORM iterations required for any problem is two. It should also be noted that since the initialization procedure to find a good starting point for FORM can take up to 2n+n performance function evaluations, Equation (25) will usually slightly under-estimate the actual number of performance function evaluations required by FORM. Preliminary testing versus MCS has shown that essentially all of the computational time required by FORM is consumed by the evaluations of the performance function. Thus, when comparing the computational efficiency of FORM and MCS, as implemented in RELAN, it is reasonable to assume that both methods have the same computational cost when the same number of performance function evaluations are undertaken. The accuracy of the MCS estimate of reliability can be determined based on the binomial distribution of the reliability estimate (a binomial proportion) and the Central Limit Theorem. As a result, the accuracy, or (1-&>)100% confidence limits, for any MCS estimate of reliability (a) where the number of realizations is sufficiently large (> 30), is given by : aÂ±z {a(\-a) I N s) (l)/2 0.5 (26) MC where z n is the value of the standard normal variate with cumulative probability level of coll a and N s is the number of MCS realizations used. Note that each MCS realization requires one MC evaluation of the performance function. The error term in Equation (26) can be equated to a maximum desirable error and then solved for N cs to obtain the number of MCS realizations M required. This requires a priori estimate of the expected reliability, however, a conservative estimate, of reliability would be 0.5 since this produces the maximum variance in the estimate of reliability and thus, the maximum number of required MCS realizations. 29 The vulnerability estimate using FORM, as described in Section 3.2.2, requires j FORM estimates of the probability of failure, where j is the number of discrete failure states. The same number of estimates of the probability of failure are also required for MCS estimates of vulnerability. However, each MCS realization can be used to evaluate all the performance functions corresponding to each of the j resistance levels. Therefore, for the same computational efficiency as a FORM estimate of vulnerability, j*N MCS realizations can be used to evaluate F the vulnerability. In other words, each MCS performance function evaluation is used to define the entire CDF of the system load, whereas each FORM performance function evaluation can only be used for defining a single point on the CDF. The accuracy of each probability of failure estimate is the same as that described for reliability. Thus, the accuracy of the FORM estimate of vulnerability cannot be independently assessed whereas the accuracy of the MCS estimate of vulnerability is dependent on the weights in Equation (7) and could be determined by the approximate joint confidence limits. However, the development of these confidence limits is beyond the scope of the present work. As discussed in Section 3.2.3, the resilience estimate using FORM requires that the probability of failure be assessed twice, for two modes of failure, in order to determine the joint probability of failure in both modes. Further, the number of performance function evaluations required in each failure mode is approximately double the number in the reliability estimation case because in the RELAN implementation of FORM all random variables are assumed to impact the performance function in each failure mode. Therefore, FORM generates two time steps worth of random variables in each failure mode and then attempts to minimize B based on these random variables even though, in this case study, half of them (one time step) are not used as inputs in each performance function evaluation. As a result, the number of random variables for the resilience estimation, considering the same sources of uncertainty as in the reliability or vulnerability estimation, is effectively n =2*n. Thus, the number of evaluations of the res performance function required for the FORM resilience estimate {N ), based on Equation (25), Fres is given by the following relationship: N Fres = 2{(2n +\) + q(2n +\)} = 2{(q+\)(4n+l)} = RES RES fo+l)(8n+2) (27) where q is the number of FORM iterations per failure mode. The same number of FORM iterations is required in both failure modes since the evaluation of the performance functions in 30 Equations (18) and (19) with the same random variables and statistics are the same problem for F O R M (except inverted) and therefore results in a pair of complement probabilities of failure. Comparing Equation (27) with Equation (25) shows that a resilience estimate using F O R M requires approximately four times as many performance function evaluations reliability estimate for the same system. Again, the accuracy of the F O R M as a single estimates for resilience can only be measured in comparison with M C S . As outlined in Section 3.2.3 the total number of performance function evaluations required for a M C S resilience estimate (N cs res) is not equal to the total number of M C S realizations used in M the first time step (NMcsresi)- Rather, it is related to the probability of failure in the first time step (pj)) by the following: NMCS res = N CS M rest + Pfl *NMCS resl = NMCS cond lpj\ + N CS M cond (28) where the number of times the second performance function (for time step two) is conditionally evaluated is NMCS cond- Thus, using NMCS cond in place of NMCS in Equation (26) yields the approximate confidence limits on the resilience value. It should be noted that for very reliable systems (i.e. low pj\) many M C S realizations will be required in the first time step to produce even a small number of NMCS cond and would therefore result in rather large confidence limits on the M C S resilience value. 3.3 GAs G A s , like other heuristic search techniques, differ from traditional optimization techniques in that they can be used for complex, discontinuous and nonlinear problems. They are unique compared to other heuristic techniques, such as simulated annealing, since they search in parallel, covering many areas of the search space by generating multiple sets of decision variables rather than a single point. Another important advantage of G A s is the fact that they generate a number of near-optimal solutions in addition to the optimal G A solution. A significant disadvantage of G A s is that their structure generally requires many more objective function evaluations than traditional optimization techniques. 31 According to Goldberg (1989), GAs are different from other more traditional optimization methods in the following ways: 1. GAs work with a coding of the decision variables, not with the variables themselves. 2. GAs search from a population of points, not from a single point. 3. GAs require only payoff (objective function) information, not trend, derivative, or other auxiliary data. 4. GAs use probabilistic transition rules, not deterministic transition rules. The third point is particularly important because this means that regardless of the complexity of the model required to generate the objective function value, as long as this value is passed to the GA, the optimization can proceed. Since the algorithm is based on genetics, the terminology used to describe the process is taken primarily from genetics. The algorithm generates an initial population of individuals at random and each individual represents a set of values for the decision variables (also called parameters in some problems). The decision variables are usually coded, most often as a series of binary digits, and the total combined set of coded decision variables are placed together to form a chromosome, string or an individual. Each individual is assigned a fitness value, which is essentially the objective function value, based on how well the individual performs. The individuals that are fit relative to other individuals in the population are then allowed to mate and exchange genetic information (pieces of their string) to produce offspring, while the other less-fit individuals are not allowed to reproduce and die off. The offspring make up the next population set or generation. This reproductive process is repeated in each generation and the population tends to evolve through the generations to contain high performing individuals. Within the general process described above there exists a myriad of simple and advanced GA techniques and GA parameter settings. The work described here will not cover all of these, rather, it will only discuss some of the more common types as well as those techniques which are available to be selected in the GA used in this study. Goldberg (1989) presents an overview of a number of the techniques not discussed here. The following GA procedures are discussed in the context of a simple GA (sGA), which generally means that the GA structure is the same as the first GAs developed by Holland (1975). Any advanced GA technique will be described as such. 32 3.3.1 Decision Variable Coding The main part of the GA usually operates on a coding of the decision variables. Binary coding is the most common coding method and is used in this work but other coding methods, such as grey (Dandy et al. 1996) and real (Cieniawski et al. 1995) coding, have been used with some success in the past. Binary coding (0 or 1 values) is best explained by way of a simple example below. A real value of a decision variable s, equal to 51, can be coded as a binary string of numbers in the following way: 1 1 0 0 1 1 -> 1*2 + 1*2 + 0*2 + 0*2 + 1*2' + 1*2Â° = 51 5 4 3 2 where the 6 digit (or bit) representation limits the number of possible values of s to 2 = 64. Due 6 to this type of coding, in a strict sense GAs are discrete optimizers. However, for all practical purposes, simply choosing a higher level of discretization approximates continuous decision variables adequately in most cases. Also, for some problems it is clearly feasible to have a discrete number of possible decision variable values that is not a power of two. In this case, the GA can be coded to select the minimum number of bits required for the desired number of possible values and check that all values of the binary coded decision variables are within their allowable ranges. For example, the GA would represent the binary coding for a decision variable with six possible values as a three-bit binary string. For each binary coded decision variable the GA checks whether the decision variable values are generated within their specified real value ranges. If the decision variable is outside of its allowable range, it is reassigned a random value within the range. A flow chart that describes the general steps of the GA is presented in Figure 3.3 and is accompanied by the following descriptions of each step. 3.3.2 Initialization A GA begins by randomly generating a starting population of individuals before entering the main generational loop. The population size is a GA parameter that should be selected to provide the population with a sufficient level of diversity but also so that computation time is not 33 too burdensome. The literature provides little in the way of guidance for the best setting of the population size and it is thus most often determined on a case by case basis. Goldberg (1989) suggests common population sizes are 30 to 100. Although GAs tend to be robust enough to find good solutions regardless of the random starting population of individuals, it is often a good idea to solve a problem a few times with different starting populations (i.e. multiple random seeds). This is especially true when trying to determine the best GA parameter settings or when trying to compare the GA results to another optimization technique since it will help to avoid lucky or unlucky results that occur because of the random initial population. Decode all individuals into decision variables GENERATIONAL LOOP REPRODUCTION Figure 3.3. Simple GA flow chart 34 3.3.3 Fitness Evaluation The next step in the GA process is to decode each binary coded solution to the corresponding real valued solution so that the fitness of each can be evaluated. GAs tend to be developed as maximizers and thus require the inverse or negative value of the objective function be used as the fitness value in minimization problems. In order to handle constraints, GAs most commonly use penalty functions to help guide the search away from infeasible solutions. A penalty function works to assign infeasible solutions an extremely poor fitness value relative to feasible solutions. Penalty functions are most often a power expression included in the fitness function in which the magnitude of the constraint violation is raised to a penalty exponent and then multiplied by a penalty multiplier. 3.3.4 Stopping Criteria After the fitness is evaluated for all population members in each generation stopping criteria are evaluated. A number of stopping criteria are possible based on convergence measures on the solutions, however the most commonly used stopping mechanism is simply a maximum generation limit. A useful and easy addition to sGAs that is incorporated into the GA used in this study is the ability to restart the GA from the last population in a previous optimization trial. The new optimization trial simply starts with the last population as the initial population instead of a randomly generated population. 3.3.5 Selection The algorithm proceeds to the reproduction phase if the stopping criterion is not met. The first step is to select the fit members of the present population who are allowed to create offspring or children for the next generation. Poor performing solutions most often do not get selected to reproduce and thus their genetic material is not carried into further generations. Goldberg (1989) outlines the biased roulette wheel, tournament selection and the ranking method as some of the possible selection techniques. Since the GA in this study uses tournament selection it is the only method discussed. In tournament selection, a pair of individuals is randomly selected from the population and the individual with the higher fitness is selected to mate with another individual 35 chosen in the same manner to create one or more offspring with a mixture of the parent's chromosomes. 3.3.6 Crossover Crossover is the process that exchanges the genetic information between two selected high fitness parents in order to produce one, or more often two, offspring. The genetic information exchanged is some portion of their binary strings. Two common types of crossover mechanisms are uniform and single-point crossover. In single-point crossover, one crossover location between any two binary bits is randomly chosen and the information to the right of this location is exchanged between parents (see Figure 3.4). Multi-point crossover simply involves randomly selecting more than one crossover location while uniform crossover implies that all potential crossover locations (locations between all binary bits) are subject to crossover. Uniform crossover allows for any possible combination of the parents genetic information. The crossover process is often random such that crossover occurs with a high probability rather than with absolute certainty. Goldberg (1989) reports common crossover probabilities ranging from 0.5 to 1.0. It should be noted that the reproduction process produces a constant population size throughout generations. Parents Offspring 1 0 1 1 0 0 1':' 111 1 1 .0 crossover site Figure 3.4. Single-point crossover 3.3.7 Mutation Mutation is the final reproduction operator applied in a simple GA. Mutation functions to assist the population from becoming trapped in a local optimum. The most common mutation method is a bitwise or jump mutation operator that, usually with a very low probability, can flip the value of any bit in any of the offspring produced by the crossover step. Jump mutation can make very significant changes to the values of the decision variables. For this reason, a newer mutation 36 operator called adjacency or creep mutation has also been used and shown to help improve optimization results (see Dandy et al. 1996; Yang et al. 1998). Creep mutation can only mutate decision variables, usually with a low probability, to an adjacent real value rather than any potential value. Both mutation operators can be used simultaneously. 3.3.8 Other Advanced GA Operators Elitism. In the elitist strategy, the best individual is brought forward without disruption from crossover or mutation. Although this operator is not necessary for the functioning of a GA, it is found useful in preventing the loss of important chromosomes. When elitism is not used, the best fitness value in a given generation can be lower than best fitness in the previous generation. If elitism is enabled the best fitness value never decreases during the evolution of the population as the best individual is always preserved. The use of elitism has become almost standard in the most recent applications of GAs in water resources (McKinney and Lin 1994; Oliveira and Loucks 1997; Mulligan and Brown 1998; Wardlaw and Sharif 1999) and is adopted in the present case study. Niching. In the basic selection processes described above there exists a possibility of mating between individuals that, although both are fit, they are very different in decision space. Combining these very different decision variables can lead to poor fitness values (Goldberg 1989). Niching, or fitness sharing, attempts to avoid mating between distant solutions. It forces individuals to share available resources by dividing the population into subpopulations of similar individuals. The goal of niching in multi-modal function optimization is to distribute the population over a number of different peaks in the search space, with each peak receiving a fraction of the population in proportion to the magnitude of that peak. Fitness sharing is also commonly used in multi-objective optimizations to stabilize subpopulations distributed along the Pareto surface. micro-GA. Another GA technique available with the GA used in this study is called the microGA. A few studies have found that this technique demonstrated faster convergence to the optimal region than a sGA (Krishnakumar 1989; Yang et al. 1998). The micro-GA technique is essentially an optimization trial using a series of sGAs with small population sizes (relative to normal population sizes in simple GAs) that are restarted when each sGA has converged. After 37 convergence is reached, the present elite solution is brought into the next generation and the remaining population members are randomly generated. Again, a sGA operates on this new population until it is converged and the process repeats itself until a maximum generation limit is reached. The randomly generated solutions inserted throughout generations can provide a similar level of overall genetic diversity as regularly sized populations in sGAs. The convergence criteria is based on the overall similarity of the entire population of individuals. Yang et al. (1998) classify a micro-GA and small-GA as those with population sizes of 5 and 15, respectively. They also modified the micro-GA technique used in Krishnakumar (1989) to include creep mutation and uniform crossover and dubbed this GA the securGA. As mentioned earlier, the best set of GA parameter values and techniques is most often assessed on a case by case basis. A small sensitivity analysis is also used in this study to determine the best GA settings and is presented in Section 5.2.1.2. 3.4 Tradeoff Surface Representation (TSR) Algorithm An ra-dimensional multi-objective solution technique, called the Tradeoff Surface Representation (TSR) Algorithm, that generates an efficient and accurate representation of convex tradeoff curves is developed and presented here. As in most other multi-objective solution techniques, the TSR Algorithm converts the multi-objective problem to a series of single objective optimization problems that must be solved iteratively. In general, each tradeoff point identified by the algorithm provides a maximize amount of information regarding the entire tradeoff surface because (1) the algorithm identifies the most representative tradeoff points and (2) the formulation of the objective function produces a maximum amount of information on the bounds of the tradeoff surface in areas of the surface between tradeoff points. Therefore, for a given number of tradeoff points, the TSR Algorithm will generate the most accurate representation of the tradeoff surface possible. Furthermore, the iterative nature of the technique is such that it may be used interactively with the decision-maker to quickly identify the preferred management solution. This is useful because multi-objective management models with large and complex environmental response models prohibit the generation of a large number of tradeoff points. 38 The TSR Algorithm is presented as follows. Section 3.4.1 presents and describes each step of the general TSR Algorithm graphically while Section 3.4.2 presents a mathematical formulation of the TSR Algorithm. Section 3.4.3 differentiates the algorithm presented here from other similar techniques and summary statements regarding the algorithm are listed in Section 3.4.4. 3.4.1 The General TSR Algorithm The general TSR Algorithm may be described by referring to the series of graphs in Figure 3.5 that show the steps in applying the algorithm to a two objective problem. The solid lines in Figure 3.5(a), (b), and (c) is the true tradeoff curve of the minimization objectives Z\ and Z and 2 the dashed lines are the approximated tradeoff curves at each step of the TSR algorithm. The algorithm can be applied to a multi-objective problem with any combination of minimization and maximization objectives. Figure 3.5(a) is a presentation of the objective space and some terminology used to describe it. The figure shows the feasible region, which is comprised of inferior (non Pareto-optimal) points and the non-inferior points making up the true convex tradeoff curve that we want to approximate by the TSR Algorithm. Also shown is the ideal point (Z* B), with co-ordinates d efined by the individual optimal values of each of the objectives A (Zi â€žâ€žâ€¢â€ž, Z miâ€ž). 2 The steps of the algorithm are as follows: Step (i). Two extreme starting points are defined (m for m-objectives). Although these can be used to bound the objective space searched by the TSR algorithm, bounds in objective space are not required by the TSR algorithm. The extreme points can be determined by minimizing each objective (Z,) individually and then calculating the value of all other objectives for the set of decision variable values that minimize objective Z, In Figure 3.5(a), points A and B are the extreme points in objective space. Step (ii). A straight line is defined by the extreme points A and B in Figure 3.5(b). The objective is then defined to find the maximum perpendicular distance from the line AB, in the direction of the ideal point (Z*AB)- In other words, the tradeoff curve is 'pushed' as far away from the straight line as possible. Point C, which lies on the tradeoff curve, is this point and the dashed line in Figure 3.5(b), line ACB, is an approximation of the tradeoff curve. 39 Step (iii). Step (ii) above is repeated between points A and C and then C and B on the tradeoff curve to produce points D and E, respectively (see Figure 3.5(c)). Note that for each additional tradeoff point generated by the TSR Algorithm after step (ii), the local ideal points, Z* AC defined by local extreme points A and C and Z* BC defined by local extreme points B and C, can be identified to aid in the mathematical formulation of the problem in Section 3.4.2. Each additional point on the tradeoff curve requires maximizing the distance from another line. The lines AC and CB are referred to as reference lines. This process can be repeated between any two adjacent tradeoff curve points until some stopping criterion is met. Step (iv). If the TSR Algorithm is used non-interactively, then the next step is to construct bounds on the location of the tradeoff curve between the tradeoff points that are identified. However, if the algorithm is used interactively, then this step would directly follow the determination of each tradeoff point. Figure 3.6(a) demonstrates the graphical determination of the tradeoff curve bounds and is described as follows. The formulation of the TSR Algorithm objective and the assumption that the tradeoff curve is convex provides the basis for which to bound the location of the entire tradeoff curve. The extreme points, A and B, provide a lower bound for each objective that is represented by the dashed lines, Z\ min and Z ,-â€ž, respectively, in Figure 3.6(a). The linearly segmented dashed line 2 m ADCEB connects the identified tradeoff points to form an approximate tradeoff curve. Notice however, that, if the problem is convex, the line ADCEB is a limit where, by the definition of convexity, all tradeoff points must either be found on or beyond, in the direction of the ideal point. The formulation of the TSR objective can be interpreted graphically as shifting the reference lines (AB, AC and CB) towards the ideal point until each shifted line is tangent to the tradeoff curve (e.g. line A C is parallel to line AC in Figure 3.6(a)). Therefore, each shifted line is a boundary beyond which, in the direction of the ideal point, no other tradeoff point can exist. Together, lines Z\ â€žâ€žâ€¢â€ž and Z2 â€žâ€žâ€¢â€ž, the shifted reference lines, and line ADCEB intersect to define the shaded area in Figure 3.6(a) where the remaining tradeoff points must be located. 41 A line Z, I mtn Shaded area bounded by dashed lines defines the possible tradeoff curve location between tradeoff points line Z 1 line A C line AC Z* AB ^ line AB' line BC' (a) A : ~^ line Z, ,â€ž,,, Lightly shaded area defines tradeoff curve location when the TSR Algorithm used while the light and dark areas combined define the tradeoff curve location when the Constraint Method is used line CD+ -Approximate tradeoff curve ( A D C E B ) line Z, ^ AB 0 line AB line BC (b) Figure 3.6. Determination of the tradeoff curve bounds using available auxiliary information in (a) step (iv) of the TSR Algorithm and (b) with the Constraint Method 42 Figure 3.6(b) is the same as Figure 3.6(a) except that the auxiliary information inherent in the Constraint Method is also shown. This information defines the boundaries of the tradeoff curve based on the Constraint Method assuming that it is used to generate tradeoff points C, D and E. For two objectives, the Constraint Method converts the problem to a single objective by making one objective a constraint and optimizing the other. In Figure 3.6(b), this would be accomplished by optimizing one objective and constraining the other to correspond to its objective values at points C, D and E. In order to generate points A and B in the Constraint Method, as with the TSR Algorithm, each objective is optimized independently. The only information inherent in the Constraint Method that can be used to bound the tradeoff curve is the lower bounds on the objectives, dashed lines Z\ â€žâ€žâ€¢â€ž and Z2 ,â€ž, and the convexity assumption. As m in Figure 3.6(a), convexity of the tradeoff curve means that the dashed line ADCEB forms a bound on the Constraint Method tradeoff curve. In addition, the convexity assumption can also be used to form additional bounds by extending the lines that connect the adjacent tradeoff points (dashed lines AD+, CD+, CE+ and BE+ in Figure 3.6(b)) until they intersect with each other or the lower bounds of the objectives. Lines AD+, CD+, CE+ and BE+ bound the tradeoff curve because the existence of another tradeoff point beyond any of them, in the direction of the ideal point, would mean the tradeoff curve is non-convex. Combining the various boundaries defines the area, shown by the sum of the lightly shaded and darkly shaded areas in Figure 3.6(b), in which the remaining tradeoff points must be located based on the Constraint Method. In Figure 3.6(b), the darkly shaded area represents the increased precision with which the TSR Algorithm identifies the location of the tradeoff curve as compared to the Constraint Method. The increased precision provides valuable information that is available at no extra computational cost since the same number of tradeoff points is generated by each method. Moreover, under the Constraint Method, the constraints are usually chosen in order to produce points on the tradeoff curve at regular intervals between the upper and lower bounds of the constrained objective. At the very least, they will be chosen at intervals such that the tradeoff points identified will not correspond to the location of the TSR Algorithm tradeoff points - unless this occurs by chance. In either case, as long as the same number of tradeoff points is generated by each technique and the problem is convex, the approximate Constraint Method tradeoff curve and corresponding bounds on the tradeoff curve will be less accurate than those generated by the TSR Algorithm. This result will be demonstrated by comparing the two methods in Section 5.2.2 for the case study multi-objective optimization results. 43 One additional procedure that follows from the generation of the tradeoff curve boundary and could be included in step (iv) of the TSR Algorithm is the development of a mean approximate tradeoff curve. Notice that in Figure 3.5(c) and Figure 3.6 the approximate tradeoff curve (ACBED) actually forms part of the boundary on the location of the tradeoff curve. Given the ability to completely bound the area in which the true tradeoff curve must lie, it is more appropriate to represent the approximate tradeoff curve with the mean values of the bounded area between any two tradeoff points as shown in Figure 3.7. Figure 3.7 is again the same example as in Figures 3.5 and 3.6 except that, as the insert shows, the mean approximate tradeoff curve is developed as the line segments DF and FC that bisect the shaded area between points D and C. In this example, when compared with the true and the approximate tradeoff curve in Figure 3.6, it can be seen that the mean approximate tradeoff is a closer representation of the true tradeoff curve. It should be clarified that the mean approximate tradeoff curve is not necessarily a more accurate representation of the true tradeoff curve than the approximate tradeoff curve, rather it is simply the average estimate of the tradeoff curve instead of the upper boundary. 44 Natural stopping criteria for the TSR Algorithm also follow from the development of the tradeoff curve bounds. Such criteria could be based on the achievement of a minimum target accuracy level for the approximate tradeoff curve. For example, one stopping criterion could be defined as follows: For point C, which is the tradeoff point of maximum distance from the reference line defined by tradeoff points A and B, IF the distance along the axis of any objective from the line AB to the point C is less then or equal to the target accuracy level, THEN one would stop generating tradeoff points between points A and C. In other words, if the stopping criterion is satisfied, we know that, for the selected objective, the maximum deviation of the true tradeoff curve from the line segments ACB is less than the distance along the corresponding objective's axis from the line AB to the point C. Other possible criteria could be based on the mean approximate tradeoff or based on minimizing the area that defines the tradeoff curve location. Of course, a maximum number of tradeoff points can also be specified as the stopping criterion. 3.4.2 Mathematical Formulation of the General m-Objective TSR Algorithm Zitzler and Thiele (1999) describe the general multi-objective problem formulation as the following: min / max Z(v) = {Z {y), Z (y), .... Z,â€ž(y)) subject to y = (y\,y , ...,y,) e Y x 2 2 s Z(y) = (Z,(y), Z (v), 2 Z (y)) e Z m (29) s where Z(y) is a vector function that maps / decision variables into m minimization and/or maximization objectives; y is a vector of decisions; Y is the decision variable space and Z is the s s objective space. In step (i) of the TSR Algorithm, extreme starting points are found for use in the objective function formulation. In step (ii), the objective of the TSR Algorithm is to find a point that is a maximum distance from anra-dimensionalreference hyperplane that is first defined by m 45 extreme points in objective space. In subsequent iterations of the TSR Algorithm, step (ii) is repeated between adjacent tradeoff surface points that act as local extreme points and define another m-dimensional reference hyperplane. The mathematical formulation is described as follows: maximize g subject to y = (yi, y, yd e Y 2 s Z(y) = (Z,(y),Z (y),...,ZJyj) â‚¬ Z 2 (30) s where g is the perpendicular distance in objective space, in the direction of the ideal point, from the vector of objective function values, Z(y), to the reference hyperplane. The reference hyperplane is given by: r*Z= c (31) where r = (n, r , ... r,) is a vector of constant coefficients; Zâ€” (Z\, Z , 2 2 Z ) is a vector of the T m unknown variables in m-dimensional objective space and c is a constant. The formulation above can be simplified slightly because maximizing the perpendicular distance from a reference hyperplane in Euclidean space is also equivalent to maximizing the scalar distances in any of the objectives from the reference hyperplane. Also, the objective must be defined such that solutions between the reference hyperplane and ideal point are assigned positive distances and solutions on the opposite side of the hyperplane are assigned negative distances. Thus, using the objective Zâ€ž which is to be minimized, a more complete formulation for the TSR Algorithm is given below and presented graphically for a two-objective problem with Z, = Z\ in Figure 3.8: maximize d- = {Z/ - Z,{y)} I ( Z , subject to y = (y\,y ,y ) 1 t 2 m a x - Z*) eY m Z(y) = (Z (y\Z (y),...,Z (y))eZ ] 2 n (32) where d, is the normalized scalar distance in objective / from the reference hyperplane to the value of the ith objective Z,(j); Z? is the value of the z'th objective at a point that lies on the 46 hyperplane found by solving the Equation (31) for z/' with the other m-l objective values at Z(y) (e.g. equation of line AB in Figure 3.8 is solved for Z/'); Z imax is the maximum value, or upper bound of Zâ€ž used to define the problem and Z* the individual minimum of objective / between the extreme points (points A and B in Figure 3.8). Normalizing the above objective relative to the range of possible values for the ith objective between the extreme points results in <i, having a maximum value of one and is a measure of how close the solution is to the local ideal point. A value of one for d indicates that the ideal point between the extreme points has been reached. t This normalization only functions to provide the same objective function value to any given solution, irrespective of the objective selected to be z'th objective in Equation (32). Z A 2 Z = r*Z < True tradeoff 2 { +c cwvc^ Z ( y ) v 7* ^ AB â€” S. < z * z,00 z X* R h ^ 1 max 7J Figure 3.8. Graphical presentation of a two-objective TSR Algorithm mathematical formulation based on Equation (32). 3.4.3 Previous Similar Attempts at Efficient Tradeoff Curve Approximation Three methods that are similar to the TSR Algorithm are proposed by Fandel (1972), as described Eschenauer et al. (1990), Das (1999) and Cohon et al. (1979). Each approach is similar to the TSR Algorithm because both have the general objective of maximizing the perpendicular distance from a hyperplane to the tradeoff surface in the direction of the ideal point. All these approaches however, are functionally different from the TSR Algorithm and the 47 approaches by Fandel (1972) and Das (1999) do not identify a representation of the entire tradeoff curve. The approaches by Fandel (1972) and Das (1999) are described in Section 2.6 while the specific differences in their approaches in comparison with the TSR Algorithm are presented below. The multi-objective approach proposed by Cohon et al. (1979) is described and then compared with the TSR Algorithm below. The method by Fandel (1972) is an interactive approach with the decision-maker and would not work as a posteriori method for identifying a representation of the entire tradeoff surface. Also, in Fandel's approach, in each subsequent iteration after the first tradeoff surface tangent point is determined, the decision-maker specifies the bounds of the objectives between which the next tangent point will be found. In contrast, after the TSR Algorithm finds thefirsttangent tradeoff point, subsequent iterations only find additional tradeoff points between two adjacent tradeoff points that were already determined by the algorithm - not the decision-maker. Even if the TSR Algorithm is used interactively, the decision-maker would specify the existing tradeoff points between which he/she is interested in finding another tradeoff point instead of a set of bounds that do not correspond to any existing tradeoff points. The approach by Das (1999) for characterizing the 'knee' of the tradeoff surface is focused on finding only a single point on the tradeoff surface that is a maximum perpendicular distance from a reference hyperplane rather than mapping the entire tradeoff surface approximation as in the TSR Algorithm. This general objective is the same as thefirsttwo steps of the TSR Algorithm. However, because the method proposed by Das (1999) is built on his previous multi-objective work with the NBI formulation (see Das and Dennis 1998), as described in Section 2.6, the mathematical formulation for finding the knee is not the same as the TSR Algorithm. Specifically, the formulation for finding the knee in Das (1999) includes the location of the origin of the normal vector on the reference hyperplane as an additional decision variable. The TSR Algorithm formulation Section 3.4.2 shows that this additional decision variable is not necessary. Also, Das (1999) presents the approach based only on a NLP algorithm and does not suggest that it is generally applicable with any optimization technique. Lastly, it is also worth noting that the NBI approach proposed by Das and Dennis (1998), which maps out an entire representative tradeoff surface, is not as accurate or efficient as the TSR Algorithm if the same number of tradeoff points are generated by each method. The reason for this is that the NBI approach strives to produce an even spread of tradeoff points while the TSR Algorithm strives to 48 iteratively find the tradeoff points that are the furthest from a hyperplane defined by surrounding tradeoff points. In addition to all of the above differences, the TSR Algorithm is distinctfromthe multi-objective approaches proposed in Fandel (1972), Das and Dennis (1998) and Das (1999) since it utilizes the auxiliary information available from successive iterations of the TSR Algorithm to first determine the bounds between which the tradeoff surface must lie and then estimate a mean approximate tradeoff surface between tradeoff points. Cohon et al. (1979) propose the Noninferior Set Estimation (NISE) Method for efficiently identifying a convex two-objective tradeoff curve in LP problems. The goal of the NISE Method is to efficiently approximate the tradeoff region in objective space within a prespecified error bound. First, the NISE Method maximizes the distance from a reference line, defined by adjacent tradeoff points, by determining the tradeoff point that is tangent to the shifted reference line. This tangent point is found by solving an LP with a weighted objective function. Cohon et al. (1979) show that the weights required to find this point on the tradeoff curve are equal to the components of the slope of the reference line. Then, the bounds of the tradeoff curve between tradeoff points are generated in the same manner as presented for the TSR Algorithm. Cohon et al. (1979) state that the NISE Method requires that the slope of the reference line is negative. This is always true for convex two-objective problems in which both objectives are to be maximized or minimized. Cohon et al. (1979) also state that the NISE Method is limited to problems with two objectives. Later work by Solanki et al. (1993) present the reasons why the NISE Method is limited to problems with only two objectives and develop a modified the NISE Method to be applicable to LP problems with more than two objectives. Although the TSR Algorithm is very similar to the NISE Method, there are four differences. First, because the TSR Algorithm was developed independent of the NISE Method, the formulations in each method for finding the maximum distance from the reference line are different. Second, although the bounds on the tradeoff curve location are determined in the same manner, the TSR Algorithm suggests the construction of a mean tradeoff curve approximation. Third, the TSR Algorithm is presented as a general method rather than being associated with one type of optimization technique as in the NISE Method. Finally, the NISE Method requires significant modifications to be applicable in problems with more than two objectives while the 49 formulation of the TSR Algorithm may not require significant modifications to be applicable to problems with more than two objectives. To the author's best knowledge, the complete TSR Algorithm as presented here is a new general approach to multi-objective optimization. However, given some of the similarities to the approaches by Fandel (see Eschenauer et al. (1990)), Das (1999), Cohon et al. (1979) and Solanki et al. (1993) it is possible that the TSR Algorithm has been previously suggested as a general multi-objective technique. In fact, if the TSR Algorithm was shown to require the same modifications as the NISE Method for application in problems with more than two objectives, then the formulations in the TSR and NISE Methods would be equivalent. The only differences would be that the TSR algorithm suggests the construction of a mean tradeoff surface and the algorithm is presented a general multi-objective technique. Therefore, more research needs to be done regarding the application of the TSR Algorithm to problems with more than two objectives and a further review of the literature is necessary to determine if the NISE Method or its variant for problems with more than two objectives have been generalized. At the very least, this work serves to re-introduce a powerful multi-objective technique to the present body of literature concerning multi-objective optimization. 3.4.4 â€¢ Summary Notes on the TSR Algorithm Although not mathematically proven, the graphical representations of the TSR Algorithm show it to be the most accurate and efficient way possible to approximate the tradeoff surface for a given number of tradeoff surface points greater than two. For example, in mathematical terms, one iteration of the TSR Algorithm, in two dimensions, applied to a convex curve between points A and B will produce point C and result in a linearly segmented approximation, line ACB, of the curve. If only one point is used between points A and B to approximate the convex curve, then point C results in the minimum area between the convex curve and any linearly segmented approximation of the curve. â€¢ The TSR Algorithm is applicable tora-dimensions.In the general w-dimensional case, the distance from a hyperplane is maximized and m points are required to define each hyperplane. The boundary of the tradeoff surface determined between tradeoff points can be constructed in a similar way as in the two objective case except that a series of intersecting hyperplanes will make up the boundary and the resulting mean approximate tradeoff will be a 50 number of linearly segmented hyperplanes. Further research needs to be conducted to determine if small modifications to the general algorithm presented here are required for application to problems with more than two objectives. The definition of the objective function at each iteration of the TSR Algorithm guarantees that the point in objective space, corresponding to the maximum perpendicular distance from a reference hyperplane in the feasible region, will lie on the tradeoff surface. This of course can only be guaranteed if the optimization technique used can guarantee convergence to the global maximum. The algorithm is not the same as minimizing the distance to the ideal point (a common goal in goal programming). The starting points for the algorithm do not have to be defined by each of the individual objective functions' extreme values. They can just as easily be defined by the tradeoff points corresponding to the minimum limits of the objectives as defined by the decision-makers. If the extreme points (objective space bounds) are chosen as inferior points then the first resulting solution generated by the TSR Algorithm will be Pareto-optimal and the most efficient accurate point relative to the defined reference hyperplane. Although not presented as such, the TSR Algorithm also lends itself to being a very useful and understandable interactive process. For example, before proceeding with each iteration in the algorithm, the decision-maker could select the regions of objective space where the algorithm is to be applied by identifying adjacent points on the tradeoff surface that bound the region in which they are most interested. This would quickly identify the preferred regions of the tradeoff surface with a minimum amount of computation. Even though the TSR Algorithm does not work for non-convex portions of the tradeoff surface the method is still useful in non-convex problems. This is because the algorithm will converge to one of the points used to define the reference hyperplane, rather than a new tradeoff point, if the tradeoff surface is non-convex in between the extreme points. Convergence to one of these points would then suggest that another multi-objective method be used in that region of the tradeoff surface such as the Constraint Method. Further, some optimization techniques may converge to solutions that lie in convex regions of the tradeoff surface even if the tradeoff surface between the extreme points is only partially convex. 51 4 WATER QUALITY MANAGEMENT IN THE WILLAMETTE RIVER BASIN The Willamette River Basin case study is used to demonstrate the novel techniques and approaches developed in this work. Section 4.1 introduces the general details regarding the case study and the corresponding water quality response model that is used. The first stage of the research consists of the estimation of the risk-based performance indicators of reliability, vulnerability and resilience with respect to meeting a water quality standard at a specified location using FORM and MCS. The case study details relevant to this stage, such as the sources of uncertainty, assumptions and method are presented in Section 4.2. The second stage develops GA-based multi-objective water quality management models for the Willamette River Basin that incorporate the performance indicators. The objectives include minimizing total basin CBOD waste treatment costs and either maximizing system reliability or resilience, or minimizing system vulnerability. The additional case study details required for this stage, such as the waste treatment cost data generation and the various water quality management model formulations, are presented in Section 4.3. 4.1 Willamette River Basin and Water Quality System Response Model The Willamette River Basin is located in northwestern Oregon and includes the state's three largest cities, Portland, Salem and Eugene (see Figure 4.1). The mainstem of the Willamette River is 300 km long, and its flow is regulated by a number of storage and re-regulation reservoirs (Leland et al. 1997). Therivermay be divided into four distinct regions, based on their hydraulic and physical characteristics (Tetra Tech 1995b). Reach I, which extends from the mouth of the river to the Willamette Falls (River Kilometre, RK, 0-42) is influenced by tides and the Columbia River; Reach II, which extends from the Willamette Falls to above Newberg (RK 42-96) is deep and slow-moving; Reach III, which extends from above Newberg to Corvallis (RK 96-208) is shallow and fast-moving, and Reach IV, which consists of the portion of the river upstream of Corvallis (RK 208-300), is also shallow and fast moving (Tetra Tech 1995b). 52 Water pollution has been an issue in the Willamette River for a number of decades. Prior to the introduction of secondary treatment requirements in the 1970s, the river experienced severe water quality problems as a result of the discharge of oxygen demanding substances from municipal and industrial point sources (Tetra Tech 1993). Since that time, there have been substantial improvements in water quality, and the current health of the river is marginal to good (Leland et al. 1997). However, pressure on water quality in the Willamette River is likely to increase in the future, as the Willamette basin is the fastest growing and most economically developed region of Oregon (Leland et al. 1997). Consequently, a DO model has been developed to help managers at the ODEQ evaluate the potential deterioration of the water quality in the river due to increased waste discharges (Tetra Tech 1995b). The mainstem of the Willamette River receives CBOD effluent from approximately 51 point sources (Tetra Tech 1995a). At present, water quality standards are defined in terms of a minimum DO level that must not be violated during critical environmental conditions (e.g. during the minimum 7-day average flow event that occurs once every 10 years) (ODEQ 1995). However, the choice of an appropriate level of protection is difficult, as there is a continuum of risk that is not well defined. For example, at DO levels between saturation and 3 mg/L, salmonids experience chronic effects of varying severity, including reductions in swimming speed, growth rate and food conversion efficiency, whereas at DO levels below 3 mg/L stream conditions are generally acutely lethal (ODEQ 1995). In addition, impacts are more severe if exposure to low levels of DO occurs more frequently and for longer periods of time (ODEQ 1995). Consequently, the analysis of risk-based performance indicators such as reliability, vulnerability and resilience may provide greater insight into the actual impact of varying DO regimes on the aquatic environment. 4.1.1 Water Quality System Response Model A QUAL2EU (Brown and Barnwell 1987) water quality response model (version 3.22) for the Willamette River was developed by Tetra Tech (1993; 1995a). This model is incorporated here as the DO response model for the Willamette River. Since the uncertainty analysis capabilities of QUAL2EU are not utilized in this study, the water quality model is only referred to as QUAL2E. The model is one-dimensional, steady-state and includes sediment oxygen demand (SOD) and average daily phytoplankton growth effects on DO. It consists of 141 model segments, each of 54 which is subdivided into computational elements of 0.16 km in length, and incorporates inflows from 14 tributaries and 51 point sources that discharge to the river. The model is calibrated using data from August 1992, verified using data from August 1994 and is considered to be valid for the summer low flow season (July to September), which is the critical period for DO. Some of the limitations of the model are that it does not incorporate the effect of periphyton production on DO, it does not account for tidal mixing with the Columbia River and it does not consider non-point or diffuse sources of nutrients or oxygen demanding substances. Furthermore, Tetra Tech (1995a) does not consider the calibrated model accurate below RK 10. QUAL2E does model nutrients such as nitrogen and phosphorous and is capable of modelling other constituents. However, since water quality problems have typically concerned low DO levels, and thus river management has been focused on mainly DO and CBOD, DO is the only constituent output analyzed from the QUAL2E model of the Willamette River. 4.2 Estimation of Reliability, Vulnerability and Resilience This portion of the case study is focused on developing the FORM and MCS based estimates of the three performance indicators and then comparing them for a range of uniform CBOD wasteload removal levels for a subset of the most important point sources. 4.2.1 System Uncertainty Before identifying the sources of uncertainty, the context and application in which the risk-based system performance indicators are going to be used need to be clearly defined. In this case study, a model is developed by combining QUAL2E with FORM and MCS to analyze the performance indicators for two scenarios. In one scenario, the focus is on the evaluation of the water quality system reliability and vulnerability with respect to meeting a DO standard at the end of the river (RK 0) during the annual low flow event. In the other scenario, the system resilience is evaluated with respect to the same DO standard at RK 0 but rather than a single event, this measure is evaluated over the entire low flow season. In each scenario the same sources of uncertainty are considered however, for some of the sources, their characterization is different. 55 As mentioned in Section 2.2, the predictions of environmental response models are subject to various types of uncertainties. The general types of uncertainty considered here are natural and parameter uncertainty. uncertainties. According to Burges and Lettenmaier (1975) these are Type II Using the classification by Tung and Hawthorn (1988), only inherent and parameter uncertainty is considered. Based on Morgan and Henrion's (1990) characterization of quantities in environmental models, the natural and parameter inputs considered uncertain here are empirical quantities. To further clarify the uncertainty in empirical quantities, according to Morgan and Henrion's (1990) classification, it is assumed here that the data used to develop the statistics for each of the uncertain quantities was obtained with no measurement error, and that the only type of uncertainty considered is variability. The variability is further assumed to be composed solely of the variability of the quantities over time and space rather than including empirical uncertainty in quantities due to incomplete scientific data. The assumptions regarding Morgan and Henrion's (1990) characterizations are clarified with respect to each source of uncertainty in the following sections. It is assumed that the QUAL2E model adequately describes all of the processes affecting DO levels in the Willamette River, although this is not strictly correct (see Section 4.1.1). As such, it is assumed that there is no significant model uncertainty. The specific sources of uncertainty considered here are the uncertainty associated with various Ka and SOD coefficients as well as the natural variation of four tributary flows and the river temperature. A description of each source and the motivation for considering them are presented in Sections 4.2.1.1-4.2.1.3. In general, a number of other sources of uncertainty such as the uncertainty associated with the hydrodynamic portion of the QUAL2E model, deoxygenation coefficient estimate and point source wasteload variability could be included in this case study. However, because the goal of this work is to demonstrate new techniques rather then use the results as the basis for actual management decisions on the Willamette River, the selected sources of uncertainty are assumed sufficient for demonstration purposes while reasonably approximating the true system uncertainty. Furthermore, when using FORM for the estimation of the performance indicators, the number of random variables included should be kept to a minimum in order to reduce computational costs. 4.2.1.1 Natural Variability: Flow and Temperature Values Flow and temperature are used to represent natural variability here, as they have been found to be important in previous studies (Takyi 1991) and there are sufficient data available to 56 characterize them. It is assumed that all measurements are exact values (i.e. no measurement errors) and that upon analyzing each flow and temperature data set, the statistical distributions and corresponding means and variances can be determined with certainty. The QUAL2E model uses inputs from the headwaters of the mainstem Willamette River and from 14 of its tributaries. On inspection of the flow records for each of these, only the headwater flows, at the confluence of the Coast and Middle Fork of the river, and the flows from the McKenzie, Santiam and Clackamas Rivers are considered as random variables, as their combined flows are approximately one order of magnitude larger than the combined flows in the remaining tributaries. The flow data used are obtained from the USGS website (www.usgs.gov) and the information related to the USGS flow measurement stations is summarized in Table 4.1. The headwater flow data are generated by adding flows at two USGS gauging stations while the McKenzie River flow is estimated based on a flow balance using the headwater flow (14157500 + 14152000) and a mainstem gauging station (14166000) just after the McKenzie River dischargers to the mainstem of the Willamette River. Each approach for flow data generation is consistent with the methods used by Tetra Tech (1993; 1995a) to determine the calibration and verification flow values. The flow data for the Santiam and Clackamas Rivers is directly available from the USGS gauging stations without further modification. Table 4.1. Information at utilized flow and temperature stations Variable USGS Station # Time Step Years of Record Flow in Willamette River at confluence of Coast 14157500+ 14152000 Daily & Middle Fork Flow in McKenzie River 1993-97 14166000- Daily (14157500+ 14152000) Flow in Santiam River 1966-86, 1988-91, 14189000 1966-86, 1988-91, 1993-97 Daily 1966-86, 1988-91, 1993-97 Flow in Clackamas River 14210000 Daily 1966-86, 1988-91, 1993-97 Temperature at Salem 14191000 Daily 1977-78, 1980-81, 1983, 1985-87 57 Although QUAL2E is capable of modeling temperature, this option is not utilized in the DO model developed by Tetra Tech (1993; 1995a). Instead, temperatures in the original model are assigned values that were measured along the mainstem of the Willamette River during the calibration and verification periods. Therefore, in order to incorporate temperature uncertainty and be as consistent as possible with the developed model, the temperatures are assumed to increase uniformly from the headwaters to the mouth of the river according to the observed relationship from the measured temperature data used for calibration and verification. In this study, this relationship is maintained while randomly varying temperature at one location. Utilizing a single location is also necessary because temperature data along the mainstem of the river are sparse. The temperature data used are obtained from the USGS and the information related to the USGS temperature measurement station is summarized in Table 4.1. Salem is chosen as the location at which temperature is varied since it is the site with the most complete record. All data analyses are carried out using only flow and temperature values from July to September, as this is the time of year for which the QUAL2E model is calibrated. For the reliability and vulnerability estimation, the seven-day moving average is obtained for all data, and the statistics for the random stream and tributary flows are obtained using the annual extreme low flow event values. The low flow event is defined as the event in which the sum of all four random streamflows is a minimum. The seven-day average temperatures that occur on the same day as this low flow event are assumed to comprise the data needed to generate the corresponding temperature statistics. The data-fitting component of the software package @RISK fwww.palisade.com) is used to determine the best-fitting distribution and the corresponding distribution parameters. The Chi-Squared test is used to analyze the fit of the distributions and the distribution with the lowest Chi-Squared test statistic is selected as the best-fitting distribution. The distribution type and its parameters, as well as the Chi-Squared test statistic for the annual low flow event streamflow, tributary flows and temperature are shown in Table 4.2. Also shown in Table 4.2 are lower bounds on the random variables that are used to ensure that the inputs to the QUAL2E model are realistic. It should be noted that RELAN automatically adjusts the PDF for each random variable so that the total probability is equal to one. Correlations between random variables are discussed later in Section 4.2.1.4. 58 ea 5 Id oo cd CU -a o c on o CO vo â€¢O X CD co u I S > 00 o o H 0 0 0 0 ON oo cd o cu -a CO CO CL) T3 cu CD X CU u o o '-4â€”* CO H OO oo t o o â€”i CO o S Â£ CD O 3 X CO 5 CD -a o o cu CD 'o o c o CU S co cd *â€”i O cd CL) fe E cd C co O Q c_r "3 CL) o c 00 ft H CL) o C cd =i cd 0 0 CN 0 0 0 0 CN CN fe O CD X cu C cd Â§ 3 S CL) CO XI -*â€”Â» o OH Cd cl In Â» ^ * -s CN o -4â€”* CO CN -a a 03 l- fe CD CD Id s CO CO cu CU co > X td'.S x cd co cd X X cd ~ cu 2 cd CU CL) M cd S co cd o o 6 u s Â« .2o o CN CU > OO & , CU C 6 co CO CN VO 00 VO vb CN CN vd CN CN â€¢Â§ X cu X cu > CU CU > O U o > CU N (D O C U O 4) O E 2 cd O o c Cu) CD > I U oo T3 E o E > s co C ^ â€¢ *-i , 1 cd O CL) -4â€”Â» o E OO td Ol H Hi U =3 cd O O 3f 3 cu c^ IH nd <D C _o cd CD -4â€”Â» ^ cd cd cd -4-J '*-X6 T-J O T3 cu > CO e M Â« tti S? fe CU X cd C cu C cd .2 x o fe T 3 CU 2 B 3 rv co SH co S cd Â« co CU M c >VH CL) CD CL) _cd 'cd > cd cu cu cd -4â€”Â» cd cd ^ CO cd cd fe B co ^ x cd CL) Id C cS CN I O o CL) o cr to rS fe CD 3 u cd |Â« cd .2 3 X !_i fe c o CU ^ cd c |g H tin co O ^ X -riS C X U-H CD <+H CL, CO CO CO CN H 59 As discussed in Section 4.2.1, unlike reliability and vulnerability, resilience is estimated over the entire low flow season and thus requires statistics that adequately characterize the tributary flows and temperature at Salem for the entire season. Since the time of travel for the Willamette River is 13 days, the data needed to characterize the flows and temperature are selected to be 13-day independent averages of daily flow and temperature (not moving averages) from July 1 to September 29 of each year. The details of this assumption are discussed further in Section 4.2.2. For each flow and temperature, every 13-day average value throughout the low flow season and throughout the years of record is combined to one data set and fitted to a single distribution. These data are analyzed and presented in the same way as the annual low flow event data and are summarized in Table 4.3. It should be noted that this treatment of the uncertain flows and temperature data assumes that the data and statistics are stationary over time. That is, the data are assumed to have no significant trend through the low flow season or years of record. Unfortunately, some of these assumptions do not hold true for all flows and temperature as shown in Appendix I. The mean and standard deviation of the flows of the Santiam River and Coast+Middle Forks tend to increase significantly in the latter half of the low flow season and therefore are non-stationary. The mean of the temperature is also non-stationary since it decreases substantially in the latter half of the low flow season. However, since the focus of this work is on demonstrating new techniques, no attempt is made to transform the data in order to remove any non-stationarities. 4.2.1.2 Parameter Variability: Sediment Oxygen Demand (SOD) Coefficient SOD is considered uncertain here because it has been found to have a marked effect on the results obtained from the QUAL2E model of the Willamette River (Tetra Tech 1995a). The SOD sampling data gathered by the modellers for calibration and verification purposes are used here to describe the SOD variability during the annual low event and over the low flow season. The use of the SOD sampling data to characterize the SOD variability in this study requires a brief presentation of the original sampling effort. 60 fl CD OJ _fl 15 > |fe ON OJ u fl OJ ^> Hâ€”Â» CJ OJ OH H a- ro CN Â© oo V O O OJ "fl VD Â© o >. C H co I OJ ON OJ CO H vo o ,fl ^ OH co 03 -3 â€ž1 OJ 2 ^ <â€” co â€¢ I-H CO g, ro 00 c o '-4â€”Â» OJ fl o c o 03 CO X) c fl o PQ CO oo oo oo t-~- CN CN CN CN â€”< CO OH OH fl fl ^Â§ OJ fl Â« o3 CD" o3 y CO co fl ~ ON 53 o3 OJ CO IT) CN Â° 2 -fl .O ID Â« CO c _o OJ fl ^ OJ ^3 O ON f Â© fl QJ lO ON o uo fl 1^ "3" OJ QJ 'â€” CD 0 0 > CD CO CO o3 C O C fl *Hâ€”Â» 03 o OH â€ž CO CN / - N CO i Q fl Â© <D OH fl CN l-H OJ Hâ€”Â» OJ fl o fl â€¢ i-H Â§3 OJ vo oo iâ€”i ro ro, ON a QJ ^H J-H Hâ€”' CX 03 OJ Hâ€”Â» OJ O O VO s 0 0 0 0 CO t> ro CN o3 >^ VO ifl ON CO ca X> o I-. o. -fl 03 (H OJ > OJ 1 12 OJ > <3 > > 03 O U o OJ CJ C OJ fl i o OJ c o o w t-< QJ Pi > OJ > s CO oj O fe OJ fl -a 2 B OJ 03 o CJ U 03 5 3 CO "S Â« 3 E ^ 03 o3 o^ -a cu CO m a IH I fl aj ^3 2 o ,"fl 2 ^ Â« (D cd C Hâ€”' +-Â» OJ CO Hâ€”Â» o 22 - CO cj â€¢Â«fl A H^ -O _o S-H . OH OJ 03 fl Â±3 x S OJ -fl OJ Â§ .2 OJ fl fl ro Â« CD CD 2 u 03 P w_ I-H =3 PH OJ O OJ c o C OJ CD OJ nfl r~* CO OJ CD OJ OH ^3 CO ^ oj >-> fl Â» o3 O CO CD s-< ^ & fl > fl O g > H-H CO CO '03 fl c o _fl fl OJ OH o o jo fe OJ 03 (-i OJ tj H â€¢3 h CN A rfl H H OJ â€¢~ CO 61 The original SOD sampling is completed over a short enough time period so that the data are assumed to be independent of time (Tetra Tech 1995a). Although the modellers also determined that the SOD data are spatially independent, the samples are divided into two groups corresponding to two river reaches, RK 0 to 42 and RK 42 to 81. The mean of each group is used in the original QUAL2E model of the Willamette to represent the SOD throughout the corresponding reach. Above RK 81, the modellers assume the SOD is negligible and set it equal to zero. Other details of the sampling effort can be found in Tetra Tech (1995a). The sampling data most likely do not adequately characterize the SOD variability over the entire low flow season due to the short sampling period. However, they are assumed reasonable here since the loss of data quality is compensated by the fact that the data are site specific. In this study, the individual samples in each group are assumed to represent a random observation of the exact (i.e. no measurement error) average reach SOD value. In accordance with the original modelling effort, the SOD is assumed deterministic and equal to zero above RK 81 and both reaches below RK 81 are characterized by separate random variables. Table 4.4 summarizes the resultant SOD statistics, lower bounds and distribution types used for each reach. 4.2.1.3 Parameter Variability: Reaeration (Ka) Coefficient Ka is included as an uncertain parameter because it has been found to be the parameter that has the most significant impact on predicted DO levels in a number of hypothetical and actual case studies (Chadderton et al. 1982; Melching and Yoon 1996). Since no site-specific measurements of Ka were obtained as part of the Willamette River QUAL2E modelling effort (Tetra Tech 1995a), the Ka uncertainty must be characterized using data available in the literature from other rivers. 62 o3 -fl CO fl o3 fl CT co cn CD 1 > H i m i VO vo Â©' Â©' VO fl â€¢ 1-H fl a" W CD Xi c o CT GO CD oj H I 1 1 CO 0 0 I O oo oq oo 0 0 CO -S 0 0 0 0 fl fl lo Â» o c o CO OH CO -4â€”Â» CD o. 3 X â€¢c >, H â€¢fl fl fl o PQ CD o cd o3 o3 o3 o3 o3 O o o Z, Â£ o Z o o Z, Z, s fl fl CO CD fl o 03 > CD Hâ€”Â» o vo CN c o Â©' <o i n < Â© IO IT) Â»â€”H CO CO Â© Â© 'â€”1 CO CO 'â€”1 Â© CO CN I â€” 1 1-H <D Â© CN Hâ€”I Hâ€”i 1â€”1 Â© Â© 1â€”I ^H Â© J ^ i ao <*>C 03 IH H^ o3 co o3 O CD Xi "fl cd 2 CD co o3 , Q O 0 0 Q O co CO nr 0 0 O CD ^ -fl G CD or I W o fl CD fc CD fc <D CD CD cd cd 03 o3 t fc to CD fl g O b cd co â€¢ ~ -O =r --3 PH -O rH +-> -H H CN < a â€¢â€” -fl fl CD t CH_ cd s g o u S , O /â€”N Â© cd CN â€¢a CO <D CD ^â€”Â» o co â€¢fl j5 CD SH â€¢ " CD fl CD id CD 1OH CD -fl Hâ€”Â» c IH CO ,o cd 03 fl fl o- 0 0 O cd cd w C+H O CO CD -o ' f l fl c3 fc -fl rfl tu td ? ^ -O <U -H "fl CD fl 03 CD O SH cd A cd co i2 fd cd JH I~Â» c2 CO <D flCL> C Iâ€” OI <-> 0 0 cd "cd I-, CD cd CN u - cd H H O .S co IH 0 0 "fl Â° -2 -S 'Â£ O CD 2 cd co fc Â£ ki CD 00 w 0 0 U Â£ -2 o? â€¢fl T3 <Â£ â€¢a 3 cd â€¢fl fl o â€¢Â£ co 13 CD 03 CD CD CD > fl -fl o CD CD CD +-> cd cd '5 H^ 03 .2 .A o IH u cd â€¢fl CD "fl -o .2 CD u CD CD -fl CO Â«J 3H CD CD ^â€”Â» Hâ€”Â» CD â€¢fl fl CO 1 J3 IS fl cd Xi CD â€¢ fl oo â€¢5 c o3 â€¢fl fl 03 . Â§g o -4â€”* -C H fl "S o cd CD -fl o I- 0 0 C o CD -O CD CD 1-H o S-H O S-H cd 'cd u C ^â€”Â» Hâ€”Â» CO 'â€¢fl Q -A CO X ! â€¢a O cd _CD c cd Hâ€”Â» CO IH ccj fl C _o td cd â€¢fl CD cd O o3 > IS u CO CO C o l s fl CD vo o CD o 'â€¢fl Xi CD CD Hâ€”Â» CD 3 "2 2 CO CD* n Â» cd H Tetra Tech (1995a) determine that, of the eight possible predictive equations for Ka in the QUAL2E model, the O'Connor and Dobbins Equation (O'Connor and Dobbins 1958) is the most appropriate for model calibration purposes. The O'Connor and Dobbins Equation predicts Ka as a function of average stream velocity and depth. The uncertainty with respect to the performance of the O'Connor and Dobbins Equation is the focus of the Ka uncertainty characterization. Melching and Flores (1999) develop a database of 371 measured Ka values and corresponding stream characteristics that is utilized to determine the variability of the O'Connor and Dobbins Equation prediction error for rivers with some of the same characteristics as the Willamette. Melching and Flores (1999) use the database to develop multiple regression equations for the prediction of Ka as a function of various measured stream characteristics. Ka is measured by a variety of gas-tracer measurement techniques and the data are collected from 166 streams in 23 US states. Ka measurements to the base e are in days" and converted to values at 20 Â°C. Stream 1 hydraulic and channel characteristics such as flow rates, velocities, depths and streambed composition are also included in the database. Melching and Flores (1999) divide the database into four groups based on flow regime (channel control or pool and riffle) and discharge (greater than 0.556 m /s and less than 0.556 m /s) for which they determine four different regression 3 3 relationships. Their attempts to compare regression predictions and Ka measurements require a logarithmic transformation (logio) of the measured and predicted Ka values because the variance in the differences between these values is not constant over all Ka values. The multiple regression equations developed by Melching and Flores (1999) are not used in this study. Instead, the database (provided by C S . Melching, personal communication) is analyzed here to develop a data set of the prediction error of the O'Connor and Dobbins Equation. The measured Ka data are assumed to be free of measurement error. For each Ka measurement in the database that has associated velocity and depth measurements, the O'Connor and Dobbins Equation is used to predict Ka values at 20 Â°C. It is observed that in many instances, the O'Connor and Dobbins Equation produces extremely poor predictions relative to the measured Ka. Therefore, the original database is reduced to include only the data for which the O'Connor and Dobbins Equation predicts Ka to within Â±60% of the measured value. Removing the data associated with these poor predictions is reasonable because Tetra Tech (1995a) determined during their calibration effort that the O'Connor and Dobbins Equation for Ka estimation best fit the observed DO data. In other words, the equation most likely performed well - or at least did 64 not produce extremely poor predictions. The +60% criterion is arbitrarily selected since no quantitative information regarding the O'Connor and Dobbins prediction error is available from Tetra Tech (1995a) and is assumed reasonable for the purposes of this case study. The reduced database is then analyzed in a manner that is consistent with Melching and Flores (1999) such that the data are divided into the same groups and the same logarithmic transformation is used to transform the measured Ka and the O'Connor and Dobbins predicted Ka. The O'Connor and Dobbins prediction error data are then generated as the difference between the transformed values of the predicted Ka and the measured Ka. Only the Ka prediction error data developed from streams with flow greater than 0.556 m /s are considered similar enough to the Willamette River. The hydraulics of the Willamette River may not be the same throughout the length of the river and are not available from Tetra Tech (1995a). Therefore, the two prediction error data groups with flow greater than 0.556 m /s, grouped by channel control or pool and riffle flow regimes, are compared by large sample estimation techniques. The results show that the means and standard deviations are not significantly different between the groups. Therefore, prediction error data of both groups are combined together to form the final data set assumed to characterize the variability in the O'Connor and Dobbins Equation prediction error when applied to streams with either a channel control or pool and riffle flow regime. The variability of the prediction error of Ka is described by four random variables that are all characterized by the same statistics derived from the final data set described in the previous paragraph. The random variables are applied to the four hydraulically distinct river reaches as described in Section 4.1. The use of four random variables for characterizing the O'Connor and Dobbins prediction error is thought to be more reasonable than assuming the prediction error would be constant over the 300 km long Willamette River mainstem. The Ka prediction error data are fitted using @RISK and the resulting statistics, lower bounds and distribution types used for each reach are presented in Table 4.4. The prediction error statistics are assumed to represent the variability of the Ka prediction error, associated with good performance of O'Connor and Dobbins Equation, during the entire low flow season and the annual low flow event. Assuming the same statistics for both scenarios may not represent reality but is reasonable given the lack of data available to proceed differently. 65 In a strict sense and in the absence of available equation performance data, the uncertainty or error in the O'Connor and Dobbins Equation would be considered model uncertainty since the equation 'models' reaeration. However, the available raw data discussed above allow for a quantitative assessment of the variability of the prediction errors associated with the equation. In this sense, the performance variability is known rather than being completely uncertain. Thus, the uncertainty associated with Ka is still thought of as variability, but more correctly, it is the variability associated with the prediction error of the O'Connor and Dobbins Equation in situations where the equation works well. 4.2.1.4 Correlations Between Random Variables Cross-correlation coefficients for reliability and vulnerability estimates are calculated from the annual low flow and temperature data while the cross- and auto- correlation coefficients for resilience are calculatedfromthe 13-day independent averages for flow and temperature during the summer low flow season. Unfortunately, the original correlation matrices based on the raw data are not positive-definite - which is a required property for producing a simultaneous set of correlated random variables. In order to develop positive-definite correlation matrices the RELAN program is modified so that the input correlations are not transformed to standard normal space according to the method of Der Kiureghian and Liu (1986). Since this modification 'corrected' only the annual low flow event correlation matrix, other measures are taken to make the resilience matrix positive-definite. The original correlations are kept unchanged as much as possible in the resilience matrix. Correlations greater in magnitude than Â± 0.7 are considered significant while all others are assumed equal to 0. After making a few other modifications to the resilience matrix, it became positive-definite. To be consistent with the approach used for the resilience matrix, only correlations of greater than about Â± 0 . 7 were included as being different than zero in the correlation matrix based on the annual low flow event. The original correlations for the annual low flow event data are presented in Table 4.5 while the modified correlations utilized to estimate reliability and vulnerability in the first stage of the research are presented in Table 4.6. The original correlations for the resilience data are presented in Table 4.7 while the modified correlations utilized to estimate resilience the first stage of the research follow in Table 4.8. 66 Table 4.5. Original correlation coefficients calculated from the annual seven-day moving average low flow and temperature data used for reliability and vulnerability estimation Variable # Variable # 1 Description 1 2 3 4 5 Flow in Willamette River at confluence of Coast and 1 -0.31 0.04 -0.01 -0.56 Middle Fork 2 Flow in McKenzie River -0.31 3 Flow in Santiam River 0.04 0.39 4 Flow in Clackamas River -0.01 0.54 0.64 5 Temperature at Salem -0.56 -0.39 -0.78 -0.60 1 0.39 0.54 -0.39 1 0.64 -0.78 1 -0.60 1 Table 4.6. Modified correlation coefficients used for the first stage reliability and vulnerability estimation Variable # Variable # 1 Description 1 2 3 4 5 1 0.0 0.0 0.0 0.0 Flow in Willamette River at confluence of Coast and Middle Fork 2 Flow in McKenzie River 0.0 1 0.0 0.0 0.0 3 Flow in Santiam River 0.0 0.0 1 0.7 -0.8 4 Flow in Clackamas River 0.0 0.0 0.7 1 -0.7 5 Temperature at Salem 0.0 0.0 -0.8 -0.7 1 67 Table 4.7. Original correlation coefficients calculated from the 13-day independent average flow and temperature data used for resilience calculation Variable # 1 Description Flow in Willamette River (t) at confluence of Coast and 1 2 3 4 5 Variable # 6 7 8 9 10 1 0.14 0.62 0.01 -0.83 0.78 0.15 0.67 0.24 -0.80 Middle Fork 2 Flow in McKenzie River (t) 0.14 3 Flow in Santiam River (t) 0.62 0.42 4 Flow in Clackamas River (t) 0.01 0.64 0.46 5 Temperature at Salem (t) -0.83 -0.26 -0.83 -0.44 6 Flow in Willamette River (t+1) 0.78 0.06 0.24 -0.19 -0.64 7 Flow in McKenzie River (t+1) 0.15 0.77 0.35 0.39 -0.11 0.14 8 Flow in Santiam River (t+1) 0.67 0.26 0.69 0.07 -0.80 0.62 0.42 9 Flow in Clackamas River (t+1) 0.24 0.60 0.59 0.74 -0.62 0.01 0.64 0.46 10 Temperature at Salem (t+1) -0.80 -0.12 -0.40 0.05 0.70 -0.83 -0.26 -0.83 -0.44 1 0.42 0.64 -0.26 0.06 0.77 0.26 0.60 -0.12 1 0.46 -0.83 0.24 0.35 0.69 0.59 -0.40 1 -0.44 -0.19 0.39 0.07 0.74 0.05 1 -0.64 -0.11 -0.80 -0.62 0.70 1 0.14 0.62 0.01 -0.83 1 0.42 0.64 -0.26 1 0.46 -0.83 1 -0.44 1 Table 4.8. Modified correlation coefficients used for first stage resilience estimation Variable # Description 1 Flow in Willamette River (t) at confluence of Coast and 5 Variable # 6 7 1 2 3 4 8 1 0.0 0.7 0.0 -0.8 0. 8 0.0 0. 7 0.0 -0.8 9 10 Middle Fork 2 Flow in McKenzie River (t) 0.0 1 0.0 0.0 0.0 0.0 0.8 0.0 0.0 0.0 3 Flow in Santiam River (t) 0.7 0.0 1 0.0 -0.8 0.7 0.0 0.7 0.0 -0.7 4 Flow in Clackamas River (t) 0.0 0.0 0.0 1 0.0 0.0 0.0 0.0 0.7 5 0.0 Temperature at Salem (t) -0.8 0.0 -0.8 0.0 1 -0.7 0.0 -0.8 0.0 0.7 6 Flow in Willamette River (t+1) 0. 8 0.0 0.7 0.0 -0.7 1 0.0 0.7 0.0 -0.8 7 Flow in McKenzie River (t+1) 0.0 0.8 0.0 0.0 0.0 0.0 1 0.0 0.0 0.0 8 Flow in Santiam River (t+1) 0.7 0.0 0.7 0.0 -0.8 0.7 0.0 1 0.0 -0.8 9 Flow in Clackamas River (t+1) 0.0 0. 0 0.0 0.7 0.0 0.0 0.0 0.0 1 0.0 10 Temperature at Salem (t+1) -0.8 0.0 -0.7 0.0 0.7 -0.8 0.0 -0.8 0.0 1 Note: t is an index for the time step. Problems with non-positive-definite matrices as described above should not generally happen when high quality data are used as the basis for the correlations. Some of the assumptions regarding the flow and temperature data may be the cause of these problems. For example, all of the random tributaries are in fact regulated by reservoirs to some degree and thus the tributary flows may not be appreciably random, especially during the annual low flow events. Further, with regard to the resilience data, it is not hard to imagine that the regulation policies of the rivers change over the low flow season. Therefore, data that have non-constant relationships throughout the low flow season may be erroneously grouped together to determine the overall correlations. The non-stationary flow and temperature data as discussed in Section 4.2.1.1 is evidence that such non-constant relationships and statistics are present in the data used. Given these problems, the most rigorous way to try and ensure positive-definite correlation matrices would be to re-analyze the assumptions regarding the randomness of the tributaries and to attempt to transform the flow and temperature data so that all of their corresponding statistics are stationary. However, since the work here is a demonstration of various techniques, it is deemed acceptable to modify the matrices directly as described above rather than attempting to analyze the data quality. In general, there seems to be no available data for characterizing the auto- or cross-correlations of the Ka or SOD values. Therefore, they are assumed independent of each other for much of the analysis. This assumption is not likely true in reality since it is reasonable to imagine that if some condition has caused the SOD coefficient and/or the Ka prediction error to be high during one time period, then the conditions are likely to persist into the next time period and thus cause the SOD or Ka prediction error to remain at a higher than normal value. The resilience is therefore also analyzed assuming arbitrary auto-correlations of 0.7 for all four Ka prediction errors and both SOD values. 4.2.2 Application of FORM to a Steady-State Response Model for Resilience Estimation The application of the basic FORM technique as presented here caters only to random variables, not random processes. As such, the technique is not well suited for dynamic problems where a time series of inputs is taken to be uncertain. The FORM technique can however, under certain circumstances, be used with a steady-state water quality response model in order to roughly 69 estimate dynamic risk-based system performance indicators such as resilience. The assumptions that this requires are discussed below. Resilience is estimated for the Willamette River over the entire low flow season. Specifically, resilience refers to the probability that, over the entire low flow period, given a set of inputs leading to system failure under steady-state conditions in one time step, the inputs in the next time step will result in the system recovering from failure under steady-state conditions. It is worth clarifying that this resilience estimate is with respect to the steady-state random variable inputs and not strictly with respect to the speed of system recovery as measured by the DO at RK 0. For example, due to the spatial differences in the model inputs, the DO at Portland may actually recover from failure before steady-state is reached. However, this can not be assessed with a steady-state model. The conditional definition of resilience used in this study does not require a time-series of inputs to be generated, rather, it only requires estimates of the lag one cross- and auto-correlations for each of the random variables. This definition allows the resilience of the system to be estimated by generating repeated two-time step random events. The water quality response model utilized here and the system being modelled limits the implementation of the resilience approach. For example, when using the steady-state DO model for the Willamette River, the resilience time step should equal the travel time of the system, as the DO model cannot respond to more rapid changes in the system. Travel time for the Willamette River, as estimated by the QUAL2E model, during the average annual seven-day moving average low flow event is approximately 13 days. Therefore, the flow and temperature statistics used for resilience estimation are based on 13-day independent averages from July 1 to September 29 of each year. The statistics for characterizing the SOD and Ka uncertainty during the annual low flow event, which were originally used for the reliability and vulnerability estimation, are assumed to also represent the statistics of their corresponding 13-day averages over the low flow season. These 13-day averages input to QUAL2E are assumed to produce approximately similar DO estimates to those that would have resulted from a dynamic estimate of the DO as a function of the 13-day time-series of inputs. Due to the length of the averaging period the assumption may not hold at all times for this 70 system but, in general, for any system, the assumption becomes increasingly more reasonable as the travel time decreases. 4.2.3 Location and Corresponding Critical Dissolved Oxygen (DO) Levels The current DO standard at RK 0 is 5 mg/L (ODEQ 1995) and is used for reliability and resilience estimation. In addition, a standard of 6 mg/L is considered. The critical DO levels at RK 0 used to bound the various failure states, the potential physical impacts of being in these failure states and the numerical indicators of severity selected for vulnerability estimation at both of the standards considered are summarized in Table 4.9. These failure states pertain to the adult life stages of cold water fish, as the reach of the river investigated is generally only used for anadromousfishpassage (ODEQ 1995). The numerical indicators of severity are assumed to be zero at DO levels above the adopted standard, as failure is defined in terms of violation of a particular DO standard. This is despite the fact that some deleterious effects can occur at higher DO levels. The set of numerical indicators of severity in Table 4.9 is chosen arbitrarily for illustration purposes, as there is no information on the quantitative impacts of being in the various failure states. As stated in Section 4.1.1, the QUAL2E model developed for the Willamette River is not considered accurate at RK 0. However, because the DO is most severe at RK 0 and since this study is only for demonstration, model predictions at RK 0 are deemed acceptable and RK 0 is used as the location for all performance indicator estimates in this study. The selection of the critical location is not necessary if MCS is used to estimate the performance indicators and may not be required if FORM is used either. For example, the minimum mainstem DO concentration could be used as the load in Equation (1) in place of the DO at a specified location. Since the conditions in this case study always lead to a minimum DO concentration at RK 0, this idea is not investigated here. 71 Table 4.9. Information about the failure states used for vulnerability estimation (ODEQ 1995) DO Standard DO Range (mg/L) (mg/L) 5,6 > standard 6 5-6 Effects Severity Index (Wi) None to minor Some avoidance 0 behavior, reduced 1 Increased avoidance behavior, reduced 2.7 cruising speed, reduced appetite 5, 6 3-5 food conversion efficiency, 5, 6 4.2.4 0-3 Acutely lethal 15.2 Wasteload Levels The QUAL2E model developed by Tetra Tech (1995a) uses the 1994 effluent wasteloads emitted by each of the 51 point sources included in the model and provides no information on any of the point source raw influent wasteloads (i.e. before waste treatment). This work is focused only on the management of CBOD since it generally has the greatest impact on DO levels. Influent CBOD wasteloads for 17 of the most important point sources is obtained from the ODEQ (Personal communication with Steve Schnurbusch and Mark Hamlin). Since these 17 point sources account for approximately 93% of the total CBOD effluent wasteload emitted by the 51 point sources in the QUAL2E model of the Willamette River, the performance indicators are estimated for a variety of uniform waste treatment levels that are applied to the 17 raw influent CBOD wasteloads. The percentage CBOD removal levels selected are 35, 50, 60, 70, 80, 90 and 95 percent of the raw influent wasteload since this range generally covers the efficiencies corresponding to primary through tertiary and chemical waste treatment of CBOD. The 17 point sources which have their CBOD effluent wasteloads modified are composed of pulp and paper mills (PPMs) and wastewater treatment plants (WWTPs) and are referred to as the managed point sources. More details on the managed point sources and the influent wasteload levels are presented Section 4.3.1. 72 4.2.5 Estimation of the Performance Indicators In this work, FORM and MCS are implemented using a modified version of RELAN (Foschi and Folz 1990), a general RELiability ANalysis software package, as described in Section 3.1.1. The performance functions for this work are defined such that ambient DO levels below the standard are evaluated to be negative (i.e. system failure). In other words, the load (Z) and resistance (R) variables in Equations (1), (10), (18) and (19) are transposed. The ambient DO levels needed to evaluate the performance function are estimated using the QUAL2E water quality response model of the Willamette River. Since the wasteloads for each managed point source and the values of the uncertain variables as generated by RELAN must be input to the QUAL2E program, the source code of QUAL2E must be modified slightly in order to be linked with RELAN. QUAL2E and RELAN are both available in the FORTRAN programming language. The QUAL2E executable model and FORTRAN source code is available on the US Environmental Protection Agency's website (www.epa.gov/ceampubl/softwdos.htm). The main changes required in the QUAL2E source code to link it to the RELAN program are: (1) reassignment of the appropriate QUAL2E variables as the random variable values generated by RELAN, (2) reassignment of the new CBOD concentrations emitted by the managed point sources in QUAL2E as a function of the desired uniform percentage of CBOD waste treatment, and (3) the transfer of the QUAL2E estimate of DO at RK 0 to RELAN. A few other changes are also made to speed up the execution time of the QUAL2E code. For example, the QUAL2E input file is only read by the overall program once, during the first performance function evaluation, rather than for every QUAL2E execution and QUAL2E execution terminates after the DO at RK 0 is calculated so that no output files are created. RELAN requires a number of parameters to be selected by the user. The FORM convergence tolerance of /?, required by the nonlinear optimization algorithm, is set equal to 0.01 unless otherwise specified, and convergence is measured in absolute terms. The maximum number of FORM iterations is set to 10 for reliability and vulnerability estimates and 20 for resilience estimates. When MCS is utilized to estimate the performance indicators, 5000 simulations are used unless otherwise specified. It should be noted that the same 5000 simulations are reproduced for all reliability and vulnerability estimations because the random seeds in RELAN are not changed. This is also true for resilience although the 5000 simulations for resilience are 73 different than the 5000 simulations for reliability and vulnerability since the statistics for the uncertain data are different. In accordance with the previously stated assumptions and methods, FORM and MCS are utilized to estimate the performance indicators for DO standards of 5 and 6 mg/L over the wasteloads specified in Section 4.2.4. The approximate 95% confidence limits for all MCS reliability and resilience estimates is presented in addition to their respective point estimates. The impact of ignoring minor correlations in the reliability matrix is assessed for both FORM and MCS by comparing reliability estimates generated using the original and modified reliability matrices presented in Table 4.5 and Table 4.6, respectively. In addition, at a standard of 5 mg/L, the sensitivity of the FORM resilience estimate to including auto-correlations in the parameters is evaluated. These additional FORM resilience results are also evaluated to determine the change in the design points over various wasteloads and are discussed with reference to the impact of relaxing the FORM B tolerance when FORM fails to converge to a /? value within the maximum number of iterations. Lastly, the efficiency and accuracy of MCS and FORM for estimation of the performance indicators will be simultaneously assessed based on the equations and discussions found in Section 3.2.4. 4.3 Multi-Objective Optimization of the Risk-Based Performance Indicators and Total Waste Treatment Cost Evaluating tradeoffs between the risk-based performance indicators and waste treatment costs using a GA-based optimization model requires the determination of waste treatment cost data, optimization formulations and appropriate GA parameters. Another relevant consideration is to minimize the computational time required to run the optimization model. All case study details, methods and assumptions previously discussed in Sections 4.1 through 4.2.3 also apply in this part of the analysis unless otherwise noted. All sources of uncertainty are characterized the same way as described previously except that a methodology for removing a number of these uncertainty sources from the analysis is developed and applied in order to improve the computational efficiency of the GA-FORM application. 74 4.3.1 Waste Treatment Cost Data In order to manage the Willamette River Basin optimally, the waste treatment costs of the point source dischargers to be managed need to be considered. This first requires identifying the characteristics of the wastes and the processes that produce them for each point source. Of the 17 managed point sources, six are PPMs and eleven are WWTPs. Table 4.10 provides the names, RK locations and if available, the influent wasteload levels or production rates for each of these facilities. The influent wasteload levels and the production rates are estimated by the ODEQ (Personal communication with Steve Schnurbusch and Mark Hamlin). Some of the influent wasteload data for the WWTPs are only available as five-day CBOD (CBOD5) values rather than five-day BOD (BOD5) values. Since these two quantities are considered equal by Tetra Tech (1993; 1995a) the same assumption is made here. For each managed point source, waste treatment costs must be associated with each level of waste treatment. Generic waste treatment costs are developed from the procedures outlined in Van Note et al. (1975) and Edward C. Jordan Co. (1979). These reports characterize the various waste treatment options along with their corresponding BOD5 removal efficiencies and total annual waste treatment costs, in US dollars, for typical municipal wastewater streams and for the pulp and paper industry, respectively. The reports only cover the waste treatment technologies available at their respective publication dates. Therefore, more recent types of waste treatment technologies developed after these dates are not considered as management options. Only the management of BOD5 is considered in this work since consistent data are unavailable for other WWTP wasteload constituents such as nitrogen and phosphorous. 4.3.1.1 Development of Wastewater Treatment Plant (WWTP) Cost Data In general, Van Note et al. (1975) estimates the total annual cost of waste treatment facilities for hundreds of combinations of waste treatment processes as a function of the influent waste stream flow rate. The report gives the annual waste treatment costs (amortized capital and operation and maintenance) for the various treatment processes, as well as the resulting percentage reductions in the waste stream constituent concentrations of BOD5, total suspended solids (TSS), and other nutrients. For each WWTP, the design influent flow rate can be used to predict cost estimates via cost equations for each treatment technology. In this case, the design flow rate is taken to be the average influent flow rate over the summer low flow season. 75 b o CD PH PH 3 O o o O O o o ND C o o r- ON o o vo o ^H OH 2 D ,9 S O H o ro Q O CQ +-Â» Q o ccu pa oa o Q o CQ (J O o .-< CN a o Pi PH CN C O _; VO NO C O NO C N C N Q VO VD CN CN CN o CQ O o o V ) OO CN ON - H t O O O CN C O C N C N O ro ro ON O CJ CN a o 00 o o X! CN i-H -H 0 0 o oo ^1- CN NO C N -3- o ON CN cu ft - O 3 PH PH H PH PH PH E-H PH PH 3 -a c u - O s Pi -H cN 3 z S cS -s u cu -o o O o -o "is H 3 O CO 4 â€¢â€”> g 'o PH 3 cs "o a. o b CU l-H PH Â« PH l-H PH PH < u > CU s cu a cu cu u 00 a u < u c o CU o U <HH o <&Â» ca CM O 3 CD z cu & o o U CU z a O o O U ^H CD PH H^ G nl PH PH 3 O -* z 'B cd I z o 3 PH .1 5 -a 3 Q o 'S cd CU o 'fe CU oo o in PH <HH o & & U â€¢4â€”1 03 00 t cu 03 oo CD O ^ w CU e )-H 3 s u (73 76 The cost data in Van Note et al. (1975) were developed for typical design concentrations of the various constituents. For BOD5, this concentration is 210 mg/L. Van Note et al. (1975) recommend redeveloping the cost equations for all waste treatment options if the actual concentration of raw influent BOD5 is significantly different than 210 mg/L. Although some of the WWTP influent BOD concentrations listed in Table 4.10 are in fact much different from this 5 value, it is beyond the scope of this work to attempt to redevelop any of the WWTP costing information. Therefore, the cost information in Van Note et al. (1975) is used without modification to determine waste treatment costs for all the WWTPs listed in Table 4.10. Van Note et al. (1975) provides equations to adjust the WWTP cost estimates to present day dollars as a function of a number of economic factors. Costs are updated here to 1999 average costs since data for all economic factors are not available at consistent times for the year 2000. The 1999 average values of the economic factors required to adjust cost estimates according to the procedure in Lence et al. (1990) are listed in Table 4.11. The service and interest factor and amortization period are assumed to be equal to the values used by Lence et al. (1990). The Federal Funds Interest rate (www,economagic.com) is assumed to be most representative of the interest rates available for municipalities to fund WWTPs. The cost of land is assumed to have a value of $5000 which is chosen to be higher than the land costs assumed by Lence et al. (1990). Labour rates for average WWTPs are assumed, as in Lence et al. (1990), to be equal to the seasonally average hourly rate (http://stats.bls.gov/blshome.html). of transportation and public utilities workers The Producer Price Index for Industrial Commodities (www.economagic.com). which was formerly known as the Wholesale Price Index when the report by Van Note et al. (1975) was originally published, is converted to have a base of 100 in 1967 to be consistent with the cost updating data and procedure used by Lence et al. (1990). The National Average Wastewater Treatment Plant Cost Index was discontinued in 1990 and therefore requires that a representative value is estimated for the 1999 average. Since construction costs are likely a major factor in the determination of this index, the trends in the Construction Cost Index are assumed to give a reasonable indication of the Wastewater Treatment Plant Cost Index. Thus, the Wastewater Treatment Plant Cost Index is extrapolated from the last recorded value in 1990 to a 1999 value based on the trend in the Construction Cost Index during the same time period. 77 Table 4.11. Summary of economic factors required to update WWTP cost estimates to 1999 values Economic Factor (1999 Avg.) Units Value years 20 Interest Rate (Federal Funds) % 4.97 Service & Interest Factor % 27.0 $/man-hr 15.68 $/acre 5000 base 1982=100 126.5 base 1967=100 395.3 base 1973=100 325.0 % 8.00 Amortization Period Labour Rate Land Cost Producer Price Index for Industrial Commodities National Average Wastewater Treatment Plant Cost Index Capital Recovery Factor Based on the economic factors listed in Table 4.11, the cost updating methods in Lence et al. (1990) and the procedures in Van Note et al. (1975), the 1999 waste treatment cost estimates for the treatment options are generated for each WWTP. Different cost estimates are generated for the WWTPs because their influent flow rates are not the same. For all WWTPs, most of the identified waste treatment options are inferior because other treatment options exist that achieves a higher level of waste treatment for equal or smaller waste treatment costs. Thus, for each WWTP, the waste treatment options are analyzed to identify the non-inferior waste treatment options that differ by at least one percent in terms of the BOD removal efficiencies. The 5 resulting number non-inferior waste treatment options for the WWTPs range from six to nine while the total annual waste treatment costs range from 0.6 to 19.2 million dollars. For each WWTP, these non-inferior solutions are used as the discrete possible values of the decision variables in the water quality management models and their costs are summarized in Appendix II. Each discrete waste treatment option has an associated BOD5 removal efficiency. In order to predict the effluent concentrations of CBOD for the selected waste treatment option (which are required as input to the QUAL2E model) the removal efficiencies are first applied to the influent BOD5 concentrations listed in Table 4.10. The resulting effluent BOD concentrations are 5 assumed equal to CBOD concentrations and then converted to ultimate CBOD concentrations 5 78 (CBOD it) using the same WWTP conversion factor (2.5) as in Tetra Tech (1993; 1995a). The u CBOD it concentrations thus become the WWTPs CBOD effluent concentrations that are used as u inputs in the QUAL2E model. The selected waste treatment option is assumed to have no impact on the various nutrient concentrations (even though this is not the case in reality) such that the effluent nutrient concentrations from each WWTP remain as defined originally in the QUAL2E model by Tetra Tech (1993; 1995a). 4.3.1.2 Development of Pulp and Paper Mill (PPM) Cost Data The unpublished report by Edward C. Jordan Co. (1979) is utilized to estimate the costs of various levels of waste treatment for each of the six PPMs considered as managed point sources. It gives total and amortized total annual average cost estimates in 1978 US dollars for various levels of improvement in wastewater quality (BOD and TSS) throughout the US pulp, paper and 5 paperboard industry. The levels of improvement in the wastewater quality are defined to be two levels of improvement with respect to in-plant production process controls, and two levels of improvement with respect to designated effluent treatment technology. The report is based on a survey of approximately 88% of the US facilities in the industry operating as of 1978. Additional sampling and analysis of many representative mills across the industry generate further information for the report. The information in the report is not updated to include more recent wastewater treatment technology options in the pulp, paper and paperboard industry. Furthermore, it is assumed here that production process technologies in the various pulp and paper mill types have not changed enough since the report was compiled by Edward C. Jordan Co. (1979) to alter the its findings significantly. In order to generate cost estimates and the resulting PPM effluent quality in terms of BOD 5 concentration, a daily production rate must be specified for each mill along with an overall classification of the mill production process/general product and various sub-classifications for other mill characteristics. The estimated daily production rates and production process/general product classifications for each of the six mills considered here are provided by the ODEQ (Personal communication with Steve Schnurbusch and Mark Hamlin) and are listed in Table 4.10. Further classification of point sources includes specifying whether the mill is an existing mill or a new point source, if an existing mill, whether the wastes are discharged directly or indirectly to the receiving waters and, in some cases, the specific product produced. All mills are categorized as existing mills and all are assumed to discharge directly to the Willamette River 79 unless the mill production process/general product classification only has mill subtypes that indirectly discharge wastes. In the cases where a specific product must also be selected, the product is chosen so that the resulting BOD5 wasteloads emitted to the Willamette River are as large as possible, so as to yield conservative estimates of the performance indicators. The three components of cost for the PPMs in Edward C. Jordan Co. (1979) are the capital, energy and operation and maintenance costs. The 1978 cost estimates for each wastewater quality improvement level across all PPM types are converted to 1999 US dollars using the economic cost indices listed in Table 4.12 and the following procedure. The amortization period and labour rates are selected to be the same as the values used for the WWTPs. The bank prime loan rate (www.economagic.com) is assumed to be most representative of the interest rate applicable to private PPMs. The CCI (www.enr.com) is used to convert 1978 capital costs to 1999 costs. The ratio of the 1999 labour rate (the seasonally average hourly rate of transportation and public utilities workers from http://stats.bls.gov/blshome.html) over the 1978 labour rate in Edward C. Jordan Co. (1979) is used to convert the 1978 operation and maintenance costs to 1999 values. The energy component costs as identified in Edward C. Jordan Co. (1979) are the electricity costs and fuel costs which are assumed to be best represented by the Industrial Sector Price for electricity and the Diesel fuel price to end users, respectively (www.economagic.com). The average ratio of the 1999 values over the 1978 values for the respective energy components is utilized to convert the 1978 estimates of energy costs to 1999 values. Table 4.12. Summary of economic factors required to update PPM cost estimates to 1999 average values Economic Factor (1999 Avg.) Units Value years 20 % 7.99 $/man-hr 15.68 Base 1913=100 6060 Electricity - Industrial Sector Price c/kWh 4.39 Diesel fuel price to end users $/barrel 24.1 % 10.18 Amortization Period Interest Rate (bank prime loan rate) Labour Rate Construction Cost Index (CCI) Capital Recovery Factor 80 Since only four options or levels of wastewater quality improvement are defined in Edward C. Jordan Co. (1979), the decision variables for each PPM have four discrete possible values, except for one PPM for which only three levels are defined. Further information regarding the physical processes involved in each of these improvement options can be found in Edward C. Jordan Co. (1979). The associated total annual amortized costs calculated for each PPM range from 0.3 to 20.5 million US dollars and are presented in Appendix II. Total annual CBOD waste treatment costs for the managed point sources range from 28.8 to 117.7 million dollars US while waste treatment costs for PPMs and WWTPs are on the same order of magnitude. For each classification of PPM the effluent flow rate and BOD concentrations are given in 5 Edward C. Jordan Co. (1979) as a function of the daily production rate. As with the WWTPs, the BOD5 concentrations have to be converted to CBOD i concentrations for input to the u t QUAL2E model. In addition, for some PPM types, the effluent flow rates vary with the selected wastewater improvement option. These effluent flow rates are changed accordingly for input to the QUAL2E model. BOD5 is again assumed equal to C B O D 5 which is then converted to CBOD uU with a PPM conversion factor (4.1) used in Tetra Tech (1993; 1995a). All other effluent constituent concentrations specified in the original QUAL2E model by Tetra Tech (1993; 1995a) are left unchanged. The generic cost curve approach applied here to generate WWTP and PPM waste treatment costs is a simplified approach. In the case of WWTPs, the types of wastewater treatment facilities at the managed point sources are to be designed as if no facilities presently exist. Similarly, for PPMs, it is also assumed that they have no pre-existing treatment facilities or process improvements. Since this is obviously not the present situation, the results of this case study should be interpreted as being relevant only to a hypothetical case study rather than to presentday water quality management in the Willamette River Basin. Advances in waste treatment technology in the WWTP and PPM industries as well as advances in PPM process technologies since the publication of the reports by Van Note et al. (1975) and Edward C. Jordan Co. (1979) are not considered in this research. Any attempt to generate more representative waste treatment cost estimates for the Willamette River Basin should focus on either more recent generic cost curve approaches or on developing the actual cost curves for each managed point source. 81 4.3.2 Computational Efficiency Improvements GAs may require a large number of fitness function evaluations before approximately optimal solutions are reached. Although GA convergence rates vary from case to case, most studies have at least on the order of one hundred GA generations and commonly on the order of several hundred to a thousand. The population size multiplied by the number of generations gives the total number of fitness function evaluations required by the GA. In this GA application, a system reliability, vulnerability or resilience estimate by FORM or MCS is required for every fitness function evaluation. Furthermore, each performance indicator estimate requires a number of evaluations of the performance function, and therefore, river DO predictions by the QUAL2E water quality response model (e.g. see Section 3.2.4). Clearly, the computational burden exacted by the combination of the GA with FORM or MCS can be very large - especially when the water quality response model is complex. Therefore, the optimization problem needs to be constructed so as to minimize the computational requirements. Three possible ways to minimize the computational requirements are to reduce the number of decision variables in the GA, to reduce the number of variables that are considered random or to simplify the water quality response model. Reducing the number of decision variables allows for convergence with fewer objective function evaluations since the size of the search space will decrease. As outlined in Section 3.2.4, the FORM computational efficiency decreases as more random variables are included in the performance function. Thus, if FORM is used as the performance indicator estimation method, reducing the number of random variables results in an approximate proportional decrease in the number of water quality response model evaluations. However, reducing the number of variables considered random in any given problem will not reduce the number of water quality response model evaluations when MCS is used to estimate the performance indicators. In this study, the QUAL2E water quality response model is not simplified or modified beyond the minor source code modifications described in Section 4.2.5. The point source information available for this work provides a natural reduction in the number of decision variables from the 51 point sources included in the QUAL2E model of the Willamette River to only 17 of the more important point sources. The number of random variables, however, is reduced by the procedure outlined below. 82 4.3.2.1 Random Variable Reduction FORM provides sensitivity information regarding the importance of each random variable with respect to causing system failure (see Equation (4)). This information is utilized here to systematically remove the least significant random variables from the analysis (i.e. make them deterministic). The general idea is to only remove the random variables from the analysis that do not significantly change the system performance indicator values obtained when all variables initially thought to be random are included as such. The impact of making certain variables deterministic must be assessed over a representative range of possible management options. In this work, the importance of each random variable with respect to the reliability estimate is assessed. The general procedure is outlined briefly below. First, with all random variables included, FORM and MCS are used to estimate reliability over a representative range of management options. This step includes all potential sources of uncertainty that are initially considered significant in the system. Then, the number of random variables is reduced by identifying the least important random variables, based on the FORM sensitivity coefficients, and assigning them deterministic values. Next, using the reduced number of random variables, the reliability estimates are again calculated with FORM over the same range of management options. Then, the new FORM-based reliability estimates are compared with the original MCS estimates (i.e. with all random variables). If these reliability estimates remain sufficiently close then this process is repeated until the reliability estimates under the reduced set of random variables are deemed unacceptably inaccurate. The original reliability correlation matrix in Table 4.5 is utilized for the random variable reduction process instead of the modified reliability matrix in Table 4.6 that is used in the first stage of the research. This is due to the fact that the reduced random variable set may produce a positive-definite reliability/vulnerability and resilience correlation matrix based on the original correlations, as calculated from the data, in Tables 4.5 and 4.7, respectively. Furthermore, when any correlated random variable is removed and made deterministic, the uncertainty data are analyzed to determine a regression relationship that is used to predict the deterministic value of the removed random variable as some function of the remaining random variables. This approach attempts to retain the mean relationships between the initial set of random variables such that the observed correlations between the removed random variables and the remaining random variables are not ignored. 83 4.3.3 Multi-Objective Water Quality Management Model Formulations The Constraint and TSR multi-objective optimization methods are compared in this work. The multiple objectives considered in this study are to minimize the total basin-wide annual cost of CBOD waste treatment, maximize the system reliability during the annual low flow event, minimize the system vulnerability during the annual low flow event, and maximize the system resilience over the entire low flow season. Due to the extreme computational burden of calculating only one of the performance indicators in an optimization framework these formulations are only developed for two objectives where one is the minimization of costs and the other is either the maximization of reliability, minimization of vulnerability or maximization of resilience. The first model is called the Reliability, Vulnerability or Resilience Constrained-Cost Minimization Model and involves minimizing the total annual cost of waste treatment while meeting a specified constraint on the selected performance indicator. The performance indicator Constrained-Cost Minimization model could include any risk-based system performance indicator and is given as follows for the case where reliability is managed: maximize \/{k (Y) + A (h (Y, Q , X s )-a ) } (33) B par s where k (Y) is the total cost of waste treatment as a function of F, the vector of discrete waste treatment options for the managed point sources in the system; h (Y, Q , X ) s par is the estimated value of reliability for a given vector of waste treatment options, water quality standard at the critical location in the stream and vector of random variable distribution parameters; Q is the s water quality standard at the critical location; X par is the vector of probability distribution parameters describing the random variables in the water quality response model; a is the s specified minimum reliability level for the system and A and B are the penalty multiplier and exponent, respectively, of the GA penalty function for the reliability constraint. The penalty is only imposed if the estimated reliability is less than the specified reliability level. The objective is inverted to become a maximization objective because the GA utilized in this work, as with most GAs, is coded as a maximization algorithm. Solving Equation (33) with a range of reliability constraints yields an approximation of the tradeoff curve. 84 The second management model is based on the TSR Algorithm as described in Section 3.4 and is called the Cost-Reliability, -Vulnerability or -Resilience TSR Model. If the model is defined with reliability or resilience as the performance indicator, as well as the minimum waste treatment cost objective, then Equation (34), based on Equation (32), describes the model considering the cost objective. maximize 4 costmax Z V^cost max â€” ^ cost j -Z (y)Y + A(Z*cost cost - Z ( )f cosl y (34) where the first term is from the TSR formulation and is the normalized scalar distance in the cost objective, d , from the reference line to the cost objective, cosl Z i(y); cos Z^ost is the value of the cost objective at a point that lies on the reference line found by solving the Equation (35) for Z * c o s t with the performance indicator objective value at Z(y); Z t max is the maximum value, or upper cos bound of Z c o s t , used to define the problem and Z* t is the minimum of the cost objective between cos the extreme points. The last term in the equation is the penalty function with A and B as defined in Equation (33) and only penalizes solutions that are outside of the extreme points. Two starting extreme points must first be identified to construct the reference line on which the first iteration of the model will be based. The extreme points are defined here by the objective space co-ordinates corresponding to the individual optimal values for the cost and performance indicator and can be represented by {Z *, Z2 ) and (Z , Z *) in objective space. These points u U cos cost 2 define the reference line: Z (y)= r*Z (y) +c cost (35) 2 where Z t(y) and Z (y) are the cost and performance indicator objectives; r is a constant cos 2 coefficient, or slope of the line, and c is a constant. The penalty function is not necessary for the solution of the Cost-Reliability, -Vulnerability or Resilience TSR Model. In this research, the penalty function is used to ensure the GA quickly explores the tradeoff curve region of interest in order to maximize computational efficiency. If for a convex problem, the penalty term in Equation (34) is removed, then solutions outside of the region defined by the reference points will not be penalized but will still have a negative 85 objective function value. The only difference may be that the GA will expend more computational effort in poor regions of the feasible space. Also, including vulnerability as the performance indicator objective in the Cost-Vulnerability TSR Model, rather than reliability or resilience, requires a minor reformulation because vulnerability is to be minimized. This can be accomplished by first defining the reference points as (Z* Z ) and (Z L cosh , Z* ), where Z is u 2 L cost 2 2 the lowest possible vulnerability value, and then modifying Equation (34) so that solutions between the reference line and the ideal point are characterized by a positive value for d . cost Solving the Cost-Performance Indicator TSR Model iteratively as described in Section 3.4.1 will yield an approximate tradeoff curve representation. Each of the above models can be used to approximate the tradeoff curve between total waste treatment costs and any one of the performance indicators. The models can be solved for any DO standard at any specified critical location or even an unspecified critical location (see Section 4.2.3). Furthermore, FORM or MCS can be used to estimate the performance indicators. The use of MCS however, requires some changes to the GA code or additional post-processing of the GA solutions as described in Section 4.3.3.1. The GA parameter values used in the solution of the water quality management models are determined by a small sensitivity analysis, described in Section 4.3.3.2, that adequately represents the characteristics of the case study. 4.3.3.1 MCS as the Performance Indicator Estimation Technique in the Water Quality Management Models The general approach taken in this study to combine the GA with MCS is to estimate the performance indicators with only a small number of MCS realizations and then exploit the GA elitism function where the best solution is carried over into the next generation without disruption. If the random sets of MCS realizations are different between generations, then, as the elite solution persists over many generations, the realizations can be continually compiled to provide an increasingly accurate estimate of the performance indicator value. This requires that a different random seed be used to generate the MCS realizations in each generation. However, each performance indicator estimate in the same generation should be based on the same set of random realizations since the fitness function values of solutions in the same generation are directly compared with one another during GA reproduction. 86 MCS is incorporated with the GA for reliability estimation in the Cost-Reliability TSR Model. A small number of MCS realizations are used for the reliability estimate (40 or 80) and the model is solved without modifying the GA source code to accumulate reliabilities of the elite solutions during the optimization. Instead, the GA output data are manually post-processed to compile reliabilities for the elite solutions throughout all generations and is then sorted to determine the optimal solution. Post-processing of the GA output is necessary because the small number of MCS realizations used generates a large amount of sampling variability in the reliability estimates that prevents the GA from continually improving on the best solutions. As a result, the elite solution in the last generation of the GA is not necessarily the optimal solution. Only the solutions with a sufficient number of compiled reliability samples are then considered further and the new fitness values are then calculated for each of these solutions based on the newly compiled and more accurate reliability estimates. The solutions are then ranked according to their fitness values to determine the optimal result for the MCS-based Cost-Reliability TSR Model. 4.3.3.2 GA Parameters The discrete waste treatment options are defined in the GA as discrete integer values that correspond to the waste treatment options available at each managed point source. The GA operates on the binary equivalents of these integers and converts them to real values as described in Section 3.3.1. A small sensitivity analysis is performed on a less computationally intensive, but still representative, version of the present case study to determine the optimal GA parameter settings for the creep mutation probability, population size, and probability of crossover when the microGA technique is and is not utilized. For the sensitivity analysis and the balance of the optimization trials in this work, the GA is set to use elitism, uniform crossover and tournament selection. For all sensitivity trials when the micro-GA technique is used, the GA program automatically resets the number of offspring to one, probability of jump mutation to zero, probability of crossover to 0.5 and turns off the niching option. For the sensitivity trials where the micro-GA technique is not utilized the number of offspring is selected as two, the niching option is used and the jump probability is set equal to the inverse of the population size. For all parameter sets in the GA sensitivity analysis three different random seeds are used to ensure the results are robust over a number of random seeds. Further, the maximum generation limit is set 87 so that a total of 5000 GA evaluations (population size x number of generations) is used across all parameter sets. The test problem used for the GA sensitivity analysis is a deterministic least-cost optimization for the present Willamette River case study. This includes the QUAL2E model of the Willamette River and the corresponding waste treatment cost data. All previously uncertain sources as described in Section 4.2.1 are set to deterministic values. The random parameter values for Ka and SOD are set equal to their mean values while the random low flow and temperature values are set equal to their independent ten-year return period values with respect to the annual low flow statistics. A modified version of Equation (34) is used as the optimization formulation since the TSR Model is to be used for the bulk of the optimization results. The objectives in the deterministic management model are to minimize the total waste treatment costs and to maximize the deterministic DO (now Z ) predicted at RK 0. No DO standard is required 2 for this formulation and the reference points for the this TSR model are determined by assessing the DO values and corresponding costs at the minimum and maximum waste treatment levels. The performance of the GA as a function of each GA parameter set is then compared in terms of the normalized distance in the cost objective between the reference line (i.e. how far the tradeoff curve is pushed out) and the optimal solution and in terms of the efficiency in reaching the optimal or near-optimal solutions. 4.3.4 Method of Evaluating the Water Quality Management Solutions In order to solve the water quality management models in Section 4.3.3 a GA is linked with the performance indicator estimation program (RELAN+QUAL2E) described in Section 4.2.5. The source code for the GA utilized here (FORTRAN GA version 1.7.1) is written by Dr. David Carroll and can be found on the Internet (www.staff.uiuc.edu/~carroll/ga.html). A waste treatment cost program is also developed that converts the various discrete options for CBOD waste treatment to their associated annual costs and effluent characteristics as described in Section 4.3.1. FORM calculates the performance indicators using the same methods and FORM parameter values as in Section 4.2.5. The sources of uncertainty considered in all management models are the reduced random variable set that is identified according to the procedure in Section 4.3.2.1. 88 Figure 4.2 is a flow chart that describes the solution approach of water quality management models. After the management objectives, performance indicator estimation technique, DO standard(s), critical location and optimization model are specified, the GA starts the optimization by generating an initial population of solutions for the first generation. Note that i and j are indexes representing the generation and population individual, respectively. A solution, yj, is selected from the population and if FORM is the performance indicator estimation technique and i > 1, the solution is checked to see if it is the elite solution from the previous generation. If it is the elite solution from generation i-1, then the evaluation of the total cost and performance indicator value are omitted and the fitness value for the solution is assigned equal to the previous fitness value of the solution. Otherwise, if yj is not the elite solution from the previous generation, FORM is not used as the performance indicator estimation technique or i=l, then the total cost and performance indicator are estimated as a function of the discrete decision variable values in solution _y. y The cost program converts the decision variable values to their corresponding waste treatment cost estimates and effluent CBOD concentration and flow rate. The CBOD effluent levels are then used by RELAN in conjunction with the random variable characterizations to estimate a single value of the performance indicator. RELAN calls QUAL2E during each evaluation of the performance function for a DO prediction at the critical location as a function of the CBOD effluent levels and the RELAN generated random variable values. After the performance indicator value is estimated, the objective function is evaluated and solution^ is assigned a fitness value according to the selected management model. The GA checks to see if the present solution is the last population member to be evaluated in generation i. If more solutions remain to be evaluated in generation i, then a new solution, y \, is selected and J+ the process described above is repeated. Otherwise, the GA checks to see if generation i is the last generation to be evaluated. If i is not the last generation, then generation i reproduces to create generation i+1 and the above process is repeated. Otherwise, the optimization model is terminated. 89 Start GA: generate random initial population, i=\,j=\ Select new solution^ from generation / yes no^ Elite solution is assigned the same fitness value as calculated in generation z'-l Cost program converts discrete treatment options to costs & CBOD effluent levels for point sources RELAN estimates the risk-based system performance indicator as a function of the new effluent levels and the random variables 4 w QUAL2E estimates DO for each evaluation of the performance with the effluent levels & random variable values Objective function is evaluated for the above calculated objectives & a fitness value is assigned to yj GA reproduces with generation i to produce generation i+\ ^2 STOP GA Figure 4.2. General optimization program flow chart The general optimization program to implement the water quality management models is utilized to generate the FORM-based cost-reliability tradeoff for the Reliability Constrained-Cost Minimization Model and the Cost-Reliability TSR Model. In addition, two MCS-Based costreliability tradeoff curves, each with a different number of MCS realizations, are generated for the Cost-Reliability TSR model. FORM-based cost-vulnerability and cost-resilience tradeoff curves are also developed from the Cost-Vulnerability TSR Model and Cost-Resilience TSR Model, respectively. All performance indicator estimates in the water quality management models are evaluated for a DO standard of 5 mg/L at RK 0. For all Cost-Performance Indicator TSR Model tradeoff curves, the cost and performance indicator objective function values required to define the starting extreme points are estimated before the optimization by setting all the managed point sources to their maximum and minimum waste treatment levels. All optimal solutions identified by the water quality management models have their nominal or approximate performance indicator estimate, as estimated by FORM or MCS, checked using 5000 MCS realizations so as to provide a more accurate representation of the true tradeoff curve. The same 5000 realizations are used to check each reliability and vulnerability estimate while each resilience estimate is also checked by a consistent set of 5000 realizations that is different from the 5000 realizations used to check the reliability and vulnerability estimates. 91 5 RESULTS AND DISCUSSION In Section 5.1 the results of the reliability, vulnerability and resilience estimation using FORM and MCS for the Willamette River case study are presented and discussed. These results form the basis for an assessment of the suitability of FORM to this particular case study for the estimation of the performance indicators relative to their estimation using MCS. An additional investigation in this section also serves to assess the importance of the Ka and SOD autocorrelations on the water quality system resilience estimates. In Section 5.2 all results pertaining to the application of the GA-based multi-objective water quality management models for the Willamette River case study are presented and discussed. Section 5.2 first presents the results of the random variable reduction process and the GA sensitivity analysis. Then, the multi-objective optimization model results required to assess the performance of the TSR Algorithm relative to the Constraint Method for multi-objective optimization are presented. The FORM and MCS comparison is then expanded to assess the relative performance of FORM and MCS within the GA-based TSR treatment cost-reliability optimization model. Finally, the tradeoff curves of treatment cost versus reliability, vulnerability and resilience, respectively, are evaluated and discussed. 5.1 Application of FORM and MCS for Performance Indicator Estimation All MCS estimates of reliability, vulnerability and resilience are made using 5000 realizations unless otherwise noted. The absolute difference between the FORM and MCS estimates of the performance indicators is the FORM error relative to MCS and is referred to as the FORM reliability, vulnerability or resilience error. All MCS estimates for reliability and resilience are bounded by their approximate 95% confidence limits as determined by Equation (26) and can be interpreted as the approximate bounds within which one is 95% confident that the true system reliability or resilience will be contained. 92 5.1.1 Accuracy of the FORM Estimates of the Performance Indicators Relative to MCS The reliability of meeting two DO standards at RK 0 is assessed over a range of uniform CBOD removal levels and is presented in Figure 5.1. Figure 5.1(a) shows the FORM and MCS estimated reliability curves for a DO standard of 5 mg/L. The FORM estimates of reliability at a standard of 5 mg/L are all lower than the MCS reliability estimates but are closest to the MCS reliabilities at higher CBOD removal levels. FORM reliability estimates are worst at lower removal levels and have a maximum FORM reliability error of 0.054 and an average FORM reliability error of 0.030. For all CBOD removal levels the FORM reliability estimates are not within the approximate 95% confidence limits for the MCS estimates of reliability with the exception of the 95% CBOD removal level. As shown in Figure 5.1(b), at a standard of 6 mg/L, FORM under-estimates reliability relative to MCS over all CBOD removal levels and FORM estimates are now closest to MCS at lower CBOD removal levels. The FORM reliability errors at a standard of 6 mg/L are slightly higher than those at a standard of 5 mg/L as the maximum FORM reliability error is 0.062 and the average FORM reliability error is 0.041. Figure 5.2 shows the FORM and MCS vulnerability curves over the same range of CBOD removal levels for DO standards of 5 and 6 mg/L. The vulnerability is measured by a unitless severity index in this case study. Therefore, the differences between the FORM and MCS vulnerability measures are difficult to assess in terms of significance. At a standard of 5 mg/L, the difference in the information provided by the reliability and vulnerability estimates is due only to the inclusion of a second failure state, with DO concentrations below 3 mg/L, in the vulnerability estimate. The observed trend in the FORM vulnerability error at a DO standard of 5mg/L, which decreases as CBOD removal levels increase, is consistent with the trend in the FORM reliability errors at a standard of 5 mg/L. The reason for this is that the probability of failing below 3 mg/L becomes extremely small at the higher CBOD removal levels and therefore has very little impact on the vulnerability estimate. At a standard of 6 mg/L, the FORM vulnerability error is more severe than the FORM vulnerability error at 5 mg/L since DO concentrations between 6 and 5 mg/L are considered as an additional failure state in the vulnerability estimate. Consequently, FORM vulnerability errors at a standard of 6mg/L are the cumulative result of three FORM probability of failure estimates instead of only two. 93 - FORM for a standard of 5 mg/L â€” A. â€¢ â€¢ â€¢ â€¢ \ - MCS for a standard of 5 mg/L â€¢ â€”â€¢â€”-A- FORM for a standard of 6 mg/L - MCS for a standard of 6 mg/L > 1.5 - > 1.0 "A. ' A . ^ \ 0.5 - 0.0 30 40 50 60 70 80 90 100 % CBOD Removal Figure 5.2. FORM and MCS vulnerability estimates at RK 0 for two DO standards Resilience is estimated using the 13-day independent average flow and temperature statistics for the entire flow season as described in Section 4.2.1.1. Resilience estimates by FORM and MCS, as well as the approximate 95% confidence limits on the true resilience value, are made over a range of CBOD removal levels for DO standards of 5 and 6 mg/L and presented in Figure 5.3. 95 0.9 30 40 50 60 70 80 90 100 , , | 80 90 100 % CBOD Removal (a) 0.3 -I 30 , 40 50 , , 60 70 % CBOD Removal (b) Figure 5.3. FORM and MCS resilience estimates at RK 0 for a DO standard of (a) 5 mg/L and (b) 6 mg/L 96 For the 5 mg/L DO standard, Figure 5.3(a) shows that while FORM resilience estimates are poor below the 80% CBOD removal level, they are much closer to the MCS predictions of resilience when removal levels are greater than or equal to 80%. Also, in contrast to reliability or vulnerability, FORM over-estimates resilience at CBOD removal levels greater than or equal to 80%). The maximum FORM resilience error is 0.136 while the average FORM resilience error is 0.070. The confidence limits for resilience are wider (i.e. less precise) than the reliability confidence limits since only a proportion of the 5000 realizations, equal to the probability of failure in the first time step, will be the actual number of realizations used to estimate resilience. In fact, the probability of failure in the first time step is so small for 90 and 95% treatment at the DO of standard 5 mg/L, less than 30 realizations are generated for the resilience estimate. Therefore, in an attempt to generate more than 30 realizations to estimate resilience, which is a requirement for the use of Equation (26), 10000 time step one realizations are used for CBOD removal levels of 90% and 95%. Only 20 failures in the first time step are observed at the 95% CBOD removal level and thus the MCS confidence limits are not plotted for this point. For over half of the CBOD removal levels at a standard of 5 mg/L, FORM fails to predict resilience values within their approximate 95% confidence limits for MCS. In comparison with the standard of 5 mg/L, the plot of resilience at a standard of 6 mg/L in Figure 5.3(b) shows the approximate 95% confidence limits for the true resilience value are greatly contracted while the FORM resilience estimates seem to be further away from the corresponding MCS estimates. The confidence limits contract because at the more stringent DO standard of 6 mg/L there are more failures in the first time step and thus more MCS realizations in the second time step with which to evaluate resilience. At the 6 mg/L standard, the maximum FORM resilience error is 0.133 and the average FORM resilience error is 0.100. 5.1.2 Complete Assessment of FORM versus MCS for Estimates of the Performance Indicators This section is included to determine how many performance function evaluations are required by MCS to estimate the performance indicators to the same level of observed FORM accuracy. Methods and equations for predicting the efficiency and accuracy of FORM and MCS for estimation of the performance indicators are presented in detail in Section 3.2.4. The efficiency of FORM for the estimation of the performance indicators, in terms of the number of performance function evaluations, becomes worse as the number of random variables in the 97 problem and FORM iterations increase (see Equations (25) and (27)). FORM accuracy is problem dependent and cannot be quantitatively predicted. The MCS efficiency and absolute accuracy have a tradeoff relationship such that if the efficiency is to be improved by using fewer MCS realizations, then the resulting level of accuracy will be decreased and vice versa (see Equation (26)). MCS accuracy is also related to the true value of the probability that is being estimated. Considering the above observations, approximate operating curves can be constructed for reliability estimation in order to visualize a comparison of the overall performance of FORM and MCS. The efficiency of FORM and MCS, as the total number of performance function evaluations, can both be plotted against the number of random variables considered. For MCS, multiple straight lines can be plotted for various levels of absolute MCS accuracy. These MCS lines are linear because the number of MCS performance function evaluations is independent of the number of random variables in the problem according to Equation (26). Figure 5.4 presents example MCS and FORM operating curves for reliability estimation where Figure 5.4(a) is for an actual reliability of 0.5 and Figure 5.4(b) is for an actual reliability of 0.9 or 0.1. It should be noted that Figure 5.4(b) can be used when the actual reliability is 0.9 or 0.1 since the approximate MCS confidence limits generated by Equation (26) are the same width for any pair of complement probabilities. The number of performance function evaluations, on a logarithmic scale, are plotted versus the number of random variables in the problem. Equation (25) is used to generate the FORM curves in both figures for two (the minimum), five and ten FORM iterations. In each figure, a 95% confidence level, the respective actual reliability values and MCS accuracy levels of 0.01, 0.03 and 0.05 are used in Equation (26) to estimate the number of MCS realizations required for the three MCS lines. By comparing Figure 5.4(a) and (b) it can be seen that for a given MCS accuracy level, the number of performance function evaluations for MCS decreases as the actual reliability deviates from 0.5, whereas the number of FORM performance function evaluations is not dependent on the reliability. The number of FORM iterations and the resultant number of performance function evaluations required by FORM are problem dependent and can only be observed. 98 - Min. FORM iterations (2) FORM for 5 iterations FORM for 10 iterations 10000 9604 MCS for accuracy within 0.01 c o 1 1067 MCS for accuracy within 0.03 TO 1000 > 0) c o 384 MCS for accuracy within 0.05 O c 3 ffi O | 100 1_ o t o 0. 10 20 40 60 80 100 120 Random variables (a) 10000 within 0.05 10 20 40 60 80 100 120 Random variables (b) Figure 5.4. Approximate FORM and MCS reliability operating curves for an actual reliability of (a) 0.5 and (b) 0.9 (or 0.1) 99 These curves can be used in a general way to assess whether FORM is more efficient than MCS for a given level of accuracy. For example, in a 60 random variable problem, if FORM converged in only two iterations and an MCS estimate of the actual reliability was found to be approximately 0.5 with a large number of realizations, then Figure 5.4(a) can be used to determine that FORM would only be more efficient than MCS if FORM'S accuracy relative to the actual reliability is within Â±0.05. Otherwise, MCS would provide a more accurate reliability estimate for the same number of performance function evaluations and would therefore be more efficient. In comparison, for the same problem and results except that the actual reliability, as estimated by a large number of MCS realizations, is found to be about 0.1 or 0.9, it can be seen in Figure 5.4(b), that FORM would only be more efficient than MCS if FORM estimated the reliability to within Â±0.03 of the true value. The method of comparing MCS and FORM by the above reliability operating curves is difficult to extend to vulnerability because of the multiple failure states involved. However, as stated earlier in Section 3.2.4, the vulnerability estimation will generally become more inefficient with respect to MCS as the number of failure states increase due to each failure state requiring one additional reliability estimation by FORM. The equations and procedures given in Section 3.2.4 for accuracy and efficiency determination of the MCS resilience estimate (Equations (26) and (28)) and the efficiency of the FORM resilience estimate (Equation (27)) can also be used to produce approximate FORM and MCS resilience operating curves. However, the accuracy of the resilience estimate by MCS is dependent, not only on the true resilience value that is being approximated, but also on the probability of failure in the first time step since this impacts the number of MCS realizations used in the conditional probability estimate. As a result, the relative performance of FORM versus MCS for resilience estimation is best evaluated using the same concepts as with the reliability operating curves but on a case by case basis instead of using general resilience operating curves. 100 The relative efficiencies of FORM and MCS are also compared for estimated reliability and resilience values in this case study. Table 5.1 summarizes the relative efficiency determination of FORM versus MCS for each reliability estimate at the 5 and 6 mg/L standards. To demonstrate the comparison, reliability operating curves in Figure 5.5 (based on Figure 5.4(a)) are used to visualize the efficiency determination for a single reliability estimate. Focusing on the reliability results for a standard of 6 mg/L and an 80% CBOD removal level, the MCS reliability is 0.547, FORM is observed to require two iterations to converge and have a FORM reliability error of 0.062. Equation (26) can be used to plot another approximate MCS operating curve (for a confidence level of 95%) at an accuracy of Â±0.062 in Figure 5.5 while the number of performance function evaluations required by FORM can be determined from the two iteration FORM operating curve at 11 random variables. In this case, FORM requires only 69 performance function evaluations while MCS requires 250 in order to estimate reliability within at least Â±0.062 of the true value in 95% of MCS trials. Thus, FORM is approximately 3.6 times more efficient than MCS. Overall, the results in Table 5.1 show that FORM is significantly more efficient than MCS for all reliability estimates in this case study. FORM is most efficient at a DO standard of 5 mg/L with an average efficiency that is approximately 16 times greater than MCS. At a DO standard of 6 mg/L, the average efficiency of FORM is reduced to be 6.7 times greater than that of MCS. The difference in relative efficiencies between the DO standards is due to the larger FORM reliability errors at the 6 mg/L standard. 101 Table 5.1. FORM and MCS relative efficiency determination for reliability DO Standard (mg/L) 35 50 60 70 80 90 95 2 2 2 2 2 2 2 Number of FORM Performance Function Evaluations (W ) 69 69 69 69 69 69 69 35 50 60 70 80 90 95 2 2 2 2 2 2 2 69 69 69 69 69 69 69 % CBOD Removal FORM Iterations MCS Reliability Estimate FORM Reliability Error (NMCS) F 5 Average - 6 - - 0.392 0.620 0.765 0.865 0.941 0.976 0.986 - 0.045 0.122 0.229 0.366 0.547 0.714 0.781 Average 1. MCS estimated using at least 5000 realizations (see Section 5.1). -m- Min. FORM iterations (2) MCS Realizations Required to Achieve Observed FORM Reliability Error 0.050 0.054 0.049 0.025 0.018 0.007 0.004 0.030 0.013 0.025 0.044 0.049 0.062 0.054 0.042 0.041 - FORM for 5 iterations 360 308 293 712 655 1899 3521 FORM Efficiency Relative to MCS ( A W AW 5.2 4.5 4.2 10.3 9.5 27.5 51.0 16.0 13.5 9.5 5.0 5.5 3.6 4.0 5.5 6.7 - 933 654 348 377 250 274 380 - - FORM for 10 iterations c o Â«-Â» ro <0 > CD C o * J u c 3 Â« U c CO EL_ o t 40 60 80 120 R a n d o m variables Figure 5.5. F O R M and MCS reliability operating curves to assess the relative efficiencies of FORM versus MCS when the actual reliability is equal to approximately 0.5 102 Repeating a similar efficiency analysis for resilience estimation yields the FORM and MCS relative efficiencies in Table 5.2. As an example, consider the CBOD removal level of 70% at a standard of 5 mg/L. Since the case study has eleven sources of uncertainty and this resilience determination requires six FORM iterations for convergence in each time step, Equation (27) predicts that at least 630 performance function evaluations are required to estimate the resilience using FORM. If the observed MCS estimates of the failure probability in the first time step (0.028) and resilience (0.887) are assumed to closely approximate the true values of each respective probability, then the total number of performance function evaluations in a MCS resilience estimate can be calculated in two steps. First, using the FORM resilience error (0.103) and the MCS resilience estimate, Equation (26) predicts that, at the 95% confidence level, approximately 37 MCS realizations in the second time step (NMcscond) are required to achieve the observed FORM accuracy level. However, due to the low probability of failure in the first time step, a much larger number of total MCS realizations are required to achieve 37 failures in time step one and therefore, 37 evaluations of the performance function in time step two. Equation (28) is then used to solve for the total number of performance function evaluations required for a MCS resilience estimate (N cs res) based on NMCS cond=37 and p/i=0.02S. In this case, an MCS M resilience estimate with at least the same level of accuracy as FORM 95% of the time requires a total of 1336 performance function evaluations. Therefore, FORM is approximately 2.1 times more efficient than MCS for resilience estimation. On average, the resilience estimation efficiency of FORM is observed to be approximately 4300 and 2 times as efficient relative to MCS for DO standards of 5 and 6 mg/L, respectively. In comparison with Table 5.1, the resilience efficiencies suggest that, with the exception of one extreme efficiency measure, FORM is more efficient relative to MCS for reliability estimation than resilience estimation. As in the reliability estimation, FORM efficiency relative to MCS is also greatest at the 5 mg/L standard. This is because of the smaller failure probabilities for the first time step and/or the smaller FORM resilience errors at the 5 mg/L DO standard. 103 --. ro â€”i oi CN - _ ro C N oo IT, ON CM ro t-m O N ^ â€”' O C O CN .3 cn 00 C N o 1-4 -fl CS â€” <D o QJ QJ -H O fl o t < r/ O to O N â€”i O N â€” ' C N N O N O N O m ro rs ID or o N O oo ro r- o o ro OO CN w â€”I -<3- ! 2 O N ^ I yI< ! 3 T NO ro CD tu i - ID > -s a ^ â€¢Â§ c 2 c3 3 o ; H g .Â§ g fc ^ I O 12 rÂ£ IT) i n o oo r o NO rN g m Tf ' o ro cn cn T-, O N in Z &2 * c (D 3 Â« S 3 pH NO ro O t C N C N CD CN ro O N O O O â€”i O â€”' â€” ro O O O O O N O O O 1^ oo OO O O O O o O CN â€”i â€”' ro ro r~- o r CN O N o O O rH Â© o CO U Di 00 3 <D _ (J rN <D ll 3 in r- NO t~r- oo r- oo C N O O O N C N r-~ o P - O O N O N O Â© O r- o â€”'< ro ro ro -rt- o o N O O fl ro N O in Â© CD +H N O O O fl s Â°O O CN -3- r-~ ro o o >0 o o CN - H o o o o o o o o o o O O O O C N o <N o o CN cv-> "tf rn â€”I ^ - M fO ^ O o r- o o CN CN ro O o o CN CN r- r- ro N O CN r-~ oo ro CN O O CN CO CD CD in irj" O H Â© O O o ^ cn to fl CD g \B â€ž ' to >H ITH UH ^ o IS cn CN CD CD oo CN O C3 Â° -s C N < N TJ- NO CN CN CN CN CN i t Et m 2 ^ g s CO O- ID tH â€” 3 Â£? C cd co ^ cS ' cd O .3 Es O m O N O O O r- oo O N O O N O N O m O N O CU Ml <u S-H CD P o o r- oo > C CD io +-â€¢ CD Â© 'fl -fl CD IH w m ro c U p to CD -*3 Q S Â° ir> Si a S Q u M <u 3 .S3 Q Some final points should be noted regarding the FORM and MCS relative efficiency determination for reliability and resilience. First, the efficiency estimates are approximations because the true reliability and resilience values are based on their 5000-realization MCS estimates and are not known exactly. Further, the resilience efficiency determination is less exact than the reliability efficiency determination because the MCS estimates of resilience have much wider approximate 95% confidence limits. Lastly, the number of performance functions required by MCS to produce equivalent FORM accuracy levels are estimated such that approximately 95% of the time MCS is used it will produce a reliability or resilience estimate that is at least as accurate or more accurate than the observed performance of FORM. This is likely to be a very conservative test in favour of FORM since 95% of the MCS trials are required to outperform or equal the FORM reliability or resilience accuracy levels. It could be argued that MCS is a more efficient method than FORM if, for the same number of performance function evaluations, MCS performed more accurately than FORM in just over 50% of MCS trials. 5.1.3 Additional Investi gations 5.1.3.1 Impact of Reliability and Vulnerability Correlation Matrix Modification As described in Section 4.2.1.4, the original correlations determined from both sets of raw flow and temperature data are not used for estimation of the performance indicators. The modifications made to the reliability/vulnerability and resilience correlation matrices are consistent with each other and required because the original resilience correlation matrix is not positive-definite. However, the original reliability/vulnerability correlation matrix is positive- definite and is therefore used to determine an additional reliability curve at a DO standard of 5 mg/L. Figure 5.6 presents the FORM and MCS reliability curves for the original and modified reliability/vulnerability correlation matrices (see Table 4.5 and Table 4.6). The difference in the FORM reliability estimates using both correlation matrices is only a small fraction of the observed FORM reliability error. Therefore, the modified correlation matrix is assumed to have a negligible impact on the reliability and vulnerability estimates. 105 0.8 FORM with raw correlations Â£> 0.7 15 .2 -m- MCS with raw correlations 1 0.6 0.5 -e- FORM with modified correlations 0.4 - B - MCS with modified correlations 0.3 30 40 50 60 70 % CBOD Removal 80 90 100 Figure 5.6. Impact of modified correlations on FORM and MCS reliability estimates for a standard of 5 mg/L 5.1.3.2 Effect of Parameter Auto-Correlation In order to investigate the importance of parameter auto-correlations, which to this point have been assumed to be zero, resilience is estimated with FORM at a DO standard of 5 mg/L across the same range of CBOD removal levels assuming all four Ka and both SOD coefficients are auto-correlated at 0.7. Figure 5.7 presents these resilience estimates along with the previous estimates of resilience at 5 mg/L for comparison. Inclusion of the parameter auto-correlations results in a very large decrease in resilience estimates over all CBOD removal levels. The minimum decrease in resilience is about 0.12 while the maximum decrease is almost 0.5. This shows the importance of determining the true magnitudes of parameter auto-correlations since even half of the observed impact would still be considered quite significant. 5.1.3.3 FORM Non-Convergence and Design Point Interpretation The additional resilience curve with parameter auto-correlations presented in the previous section also provides the results for the final additional investigation. In order to generate the final three resilience points for this curve at 80, 90 and 95% CBOD removal the /? tolerance of 106 FORM must be relaxed due to non-convergence in the first failure mode at the regular f5 tolerance setting of 0.01. The fact that FORM does not converge here is troubling since the reliabilities eventually determined at these points are not extraordinarily high or low (i.e. one or zero) with reliability values between 0.96 and 0.99. Increasing the maximum number of FORM iterations is found to be unsuccessful in achieving convergence and thus only the /3 tolerances can be changed to generate a FORM estimate of resilience. However, it is observed that if the change in /? tolerance is too great between adjacent CBOD removal levels, FORM can converge to different design points and therefore produce resilience estimates that are greatly misleading. 0.2 J 30 , 40 , 50 , 60 , 70 , 80 , 90 I 100 % CBOD Removal Figure 5.7. Resilience at a DO standard of 5 mg/L with and without parameter autocorrelations In order to demonstrate this result Figure 5.8 shows the changes in the first failure mode design point values for select random variables over the range of CBOD removal levels. The design point values are presented relative to those at the 35% CBOD removal level. Only four random variables are presented here to maintain the clarity of the figure. Relaxation of the /? tolerance to 0.05 at a removal level of 80% produces design points and a resilience estimate that are 107 consistent with the trends established at lower removal levels. However, at the 90% CBOD removal level, FORM does not converge until a f3 tolerance of at least 0.15 is selected. The variable values at the design points change dramatically in value although the actual resilience prediction at the 90%> CBOD removal level, as evidenced in Figure 5.7, continues to follow the same trend established at lower removal levels. For resilience estimation at the final CBOD removal level of 95%> an initial /? tolerance of 0.15 resulted in non-convergence and a subsequent /? tolerance of 0.25 produced the completely different design point values shown by the dashed lines and hollow data points in Figure 5.8. Only after iteratively modifying the /? tolerance, a value of 0.225 is found where FORM converges to design points that follow the similar pattern established at the lower CBOD removal levels. This is significant because the resilience estimate at a /? tolerance of 0.25 is found to decrease sharply to 0.56 rather than continue to improve, as CBOD removal levels increase, to 0.73 when the /? tolerance is 0.225. 1.50 Coast+Middle fork flow - Santiam flow Reach I Ka prediction error -0.50 Reach I SOD 40 50 - Dashed lines with hollow data points show the design point if the Beta tolerance is changed too much from the tolerances used to at the preceeding % -CBOD removal level 60 70 90 100 % CBOD Removal Figure 5.8. Normalized design points for time step 1 in resilience estimate at a DO standard of 5 mg/L with parameter auto-correlations for select random variables 108 The observed sensitivity to /? tolerance is not an attractive feature of the FORM technique since it has the potential to completely change the shape and trend in the resilience curve. In this case the convergence to a different design point seems likely to have been caused by changes to the j5 tolerance. However, this is not necessarily the only possible cause of the problem. As long as the performance function changes slightly then the possibility exists for FORM to find a drastically different design point. This type of behaviour could prove difficult to handle when the goal is to produce an incrementally changing resilience tradeoff curve. 5.1.4 Summary The results in Section 5.1.1 show that FORM predictions of the performance indicators are subject to various magnitudes of error, the performance of FORM can change between various CBOD removal levels and FORM is most accurate at a DO standard of 5 mg/L and at the higher values of reliability and resilience. Therefore, for the same standard and CBOD removal level, the FORM-based reliability estimate, for example, may be estimated with much less error than the FORM-based resilience estimate. These findings also demonstrate that FORM performance is problem dependant and serve as a reminder that changing a water quality management parameter, such as the water quality standard or the waste removal level, is equivalent to creating a new performance function and thus a new problem to be evaluated with FORM. FORM performance for estimating reliability, vulnerability and resilience for the Willamette River case study is significantly inaccurate when compared to most 5000 realization MCS estimates of the performance indicators. However, the comparison of FORM efficiencies relative to MCS in Section 5.1.2 shows that FORM is much more efficient than MCS. Given this efficiency gain, the significant inaccuracies of FORM, as well as the convergence and multiple design point identification problems, are accepted here in order to use FORM as the base technique for reliability, vulnerability and resilience estimation in the water quality management models. 109 5.2 Application of the G A to Evaluate the Multi-Objective Water Quality Management Models All tradeoff curves presented in this section are developed over the complete range of possible discrete waste treatment options for all 17 managed point sources under a DO standard of 5 mg/L at RK 0. In all cases, five points on the tradeoff curves are presented. The extreme tradeoff curve points are generated in the absence of the multi-objective water quality management models. This is achieved by selecting the minimum and maximum cost waste treatment options for all point sources and then estimating the resultant total waste treatment costs and the system reliability, vulnerability or resilience values for each extreme. The remaining three tradeoff points are generated by one of the multi-objective water quality management models described by Equations (33) and (34). 5.2.1 Preliminary Analyses for Specification of the Optimization Models 5.2.1.1 Random Variable Reduction The random variable analysis for removing the least important random variables should be carried out over a range of reliability values similar to the range that will result from all possible water quality management options. Assessment of the system reliability at RK 0 for a DO standard of 5 mg/L when all of the managed point sources treat their wastes at the maximum and then minimum CBOD waste removal level defines the range of possible reliabilities for the Willamette River. This reliability range coincides very closely with the range of reliabilities observed for the 5 mg/L DO standard in Figure 5.1(a). Therefore, the same uniform CBOD waste treatment levels used to develop the previous reliability curves are assumed to be adequately representative of the range of possible discrete waste treatment options as determined in Section 4.3.1 and are used to reduce the number of random variables considered. Also, the original reliability and vulnerability correlation matrix (Table 4.5) is used for this analysis. 110 Figure 5.9 shows the system reliability curves over a range of uniform CBOD removal levels at the end of each stage in the random variable reduction. Random variables that are shown to be the least important at each stage of the analysis are set equal to deterministic values. The importance of each random variable is assessed for each reliability estimate using the square of the direction cosines, which represent sensitivity indices for the random variables, as calculated by Equation (4). The variables with the smallest values across the majority of CBOD removal levels are successively made deterministic. In the first stage, when all eleven variables are considered uncertain, the three random variables representing the uncertainty of the Ka prediction error in the three most upstream reaches do not significantly impact the reliability estimate. Setting these equal to their mean values and estimating new reliability values for the reduced set of eight random variables produces the FORM reliability curve for the eight variable case which is nearly identical to the FORM reliability curve when eleven variables are considered uncertain. In the next stage, examination of the new set of squared direction cosines shows that the temperature random variable is the least significant. Temperature is converted to a deterministic value via a regression relationship with the sum of the four random tributary flows. The resultant FORM reliability curve for the seven variable case remains quite similar to the original eleven variable FORM reliability curve with an average absolute deviation from the original reliability curve of only 0.005. The next least significant random variable is then assessed to be the Clackamas River flow and is therefore converted to a deterministic value via a regression relationship with the sum of the remaining three random tributary flows. The resulting reliability estimates yield the FORM reliability curve for the six variable case. Again, this reliability curve is very close to the original FORM reliability curve. However, an assessment of the six groups of squared direction cosines shows that it is not possible to further identify only one random variable that is the least important over the majority of the CBOD removal levels. Therefore, the random variable reduction process is stopped with six uncertain sources remaining for inclusion in the uncertain optimization models. Ill 1.0 0.9 0.8 â€¢ M C S with 11 variables >. 0.7 .Q m .2 1 FORM with 11 variables 0.6 - FORM with 8 variables - 3 Ka are removed 0.5 - X â€” FORM with 7 variables temperature is removed 0.4 â€¢A - - FORM with 6 variables Clackamas River flow is removed 0.3 30 40 50 60 70 80 90 100 % CBOD Removal Figure 5.9. Reliability curves for various stages of the random variable reduction process The reduced set of random variables includes only the most downstream Ka prediction error, both SOD coefficients and the three tributary flows of the Middle+Coast Fork of the Willamette, McKenzie and Santiam Rivers. These six sources of uncertainty are assumed to adequately characterize the uncertainty in the Willamette River QUAL2E model and are thus used for all subsequent performance indicator estimations within the water quality management models. Also, when the original correlation matrices (see Table 4.5 and Table 4.7) are reduced to include only the six random variables, both become positive-definite matrices. Therefore, these reduced matrices (see Appendix III), with correlations originally estimated from the annual low and the 13-day independent average flow data, are utilized as the respective correlation matrices for the optimization trials. The random variables for Ka and SOD are left uncorrelated throughout the optimization trials and the initially considered random variables removed from consideration are set to the deterministic values as outlined in the previous paragraph. The motivation to reduce the number of sources of uncertainty in a problem is not an ideal characteristic of FORM since it may begin to oversimplify the case study considered. Further, in 112 some case studies it may be difficult to identify the least important random variables over the representative range of management alternatives since the least significant random variables may change considerably between alternatives. The random variable reduction process as presented here is restricted to be applicable to management situations that consider only reliability at a DO standard of 5 mg/L. Extending these results to be applicable in management situations where the other risk-based performance indicators such as resilience or vulnerability are of interest requires the assumption that the relative importance of each random variable in the estimation of vulnerability or resilience would be same as in the reliability case. Although this assumption is reasonable for the purposes of this case study, it would likely be unreasonable in an actual management situation. Lastly, this approach requires the assumption that an insignificant random variable converted to a deterministic variable would remain insignificant if it was included as a random variable in further reduced sets of random variables. The fact that MCS provides no such motivation to artificially reduce the number of random variables is a significant advantage of MCS over FORM. 5.2.1.2 GA Sensitivity Analysis A small sensitivity analysis is undertaken to determine a robust set of GA parameters to be used for the optimization trials as described in Section 4.3.3.2. The sensitivity analysis can be broken into two categories that differ in the use of the micro-GA technique and the use of regular GA parameters (i.e. no micro-GA). When the micro-GA technique is utilized, all parameter combinations of population sizes of 5 and 10 and probability of creep mutation values of 0.12 and 0.24 are evaluated for three random seeds. When only regular GA parameters are utilized, all parameter combinations of population sizes of 20 and 50, probability of crossover of 0.5 and 0.75 and probability of creep of 0.12 and 0.24 are also evaluated for three different random seeds. The results of the sensitivity analysis show that one parameter set in each main category is clearly the best of the respective categories. The best micro-GA parameter set is a population size of 5 and a probability of creep equal to 0.12 while the best regular GA parameter set is a population size of 20, probability of crossover of 0.75 and a probability of creep of 0.12. For each random seed in both parameter sets the evolution of the elite solution over the number of GA evaluations is presented in Figure 5.10. This figure clearly shows the micro-GA parameter set is the best of all of the parameter sets examined in this analysis since for all three random 113 9 seeds the global optimum for the entire sensitivity analysis is identified. In contrast, the best regular GA parameters fail to find the global optimum for any random seed. Further, the microGA parameter set produces the best fitness values across all GA evaluations. 0.58 0.56 i , 0.54 ii CD CO > m in o 0.52 c 0.5 4i One random seed using the best regular GA parameters 0.48 0.46 - One random seed using the best Micro-GA parameters r 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Number of G A Evaluations Figure 5.10. Sensitivity analysis results for each of the best micro-GA and regular GA parameter sets These results are not completely unexpected given that the micro-GA parameter set uses four times as many generations as the regular GA parameter set to develop good performing solutions. Recall from Section 4.3.3.2 that the reason for a different number of generations between parameter sets is that the maximum number of GA evaluations is set to 5000 for all GA parameter sets. As a result, the number of GA generations is a function of the population size. Although each GA technique is designed to work with different population sizes the best regular GA parameters are also tested at additional population sizes of 5 and 10. These population sizes in combination with the other regular GA parameters did not improve the results. The sensitivity analysis here must represent the conditions of the optimization trials for the real case study. The actual optimization trials will be too computationally expensive to do more than 5000 GA evaluations and thus the sensitivity analysis is only used to determine the GA parameters that produce the best results with the lowest number of GA evaluations up to a 114 maximum of 5000. For example, it is possible that if the best regular GA parameter set is used for the same number of generations as the best micro-GA parameter set, the resulting optimal objective function values may be the same or better than the micro-GA technique optimal values. However, that result would be irrelevant to the case study because at four times the allowable computational cost this parameter set would simply be too burdensome to be able to generate a complete set of results in a reasonable amount of time. The complete set of GA parameters determined by the sensitivity analysis include the use of the micro-GA technique, a creep mutation probability of 0.12, a population size of 5 and a resulting maximum number of generations of 1000 in accordance with the 5000 GA evaluation limit. The additional GA parameter values as outlined in Section 4.3.3.2 are the use of elitism, uniform crossover with a probability of 0.5, one offspring per pair of parents, a jump mutation probability equal to zero and tournament selection. The number of discrete options for each of the 17 managed point sources are available in Appendix II while the total number of possible solutions in the search space is 8.1xl0 (or 8100 billion). These GA parameters are used for the balance 12 of the optimization trials unless otherwise specified. 5.2.2 Application of the TSR and Constraint Methods for Multi-Objective Optimization Cost and reliability are used as the objectives to compare the tradeoff curves developed by the Constraint and TSR Methods. In addition, only the FORM estimates of reliability are used in this comparison rather than a MCS check on the optimal reliability since the differences between each multi-objective method are to be demonstrated independent of the reliability estimation technique. Since it is not known with certainty if any of the tradeoff curves will be convex, for all tradeoff curves, any auxiliary information presented regarding the bounds of the tradeoff curve location will not assume convexity. The non-convex uniform waste treatment level- reliability curve in Figure 5.1(b) provides a rationale for not making the convexity assumption. Thus, only the approximate tradeoff curve and partial bounds on the tradeoff curve, as developed in Section 3.4.1, are presented. 115 Figure 5.11 shows the approximate cost-reliability tradeoff curves developed by the TSR and Constraint Methods for a standard of 5 mg/L at RK 0 and the respective partial bounds on the tradeoff curve that can be constructed for each method. The endpoints of each tradeoff curve are the same and occur at the minimum and then maximum waste treatment cost options for each point source. The minimum reliability levels specified for the Constraint Method tradeoff curve points are selected so that the tradeoff points are equally spaced along the reliability axis. Three iterations of the TSR Algorithm produce the three interior tradeoff curve points for the TSR tradeoff curve. The tradeoff curves produced by each method are very similar at lower reliability levels but quite different at higher reliability levels. For example, if a decision-maker used the tradeoff curve based on the Constraint Method, they may think that approximately 100 million $US/yr would be required to reach a reliability level of 0.9. In contrast, if the decision-maker used the tradeoff curve based on the TSR Algorithm, they would see that reaching the 0.9 reliability level may instead require 80 million $US/yr. The bounds on the tradeoff curve location towards the ideal point are also shown in Figure 5.11 as dashed lines for both multi-objective methods. In the absence of the convexity assumption, the only bounds on the tradeoff curve location that can be constructed for the Constraint Method are those that use the assumption that cost and reliability are conflicting objectives. These bounds are easy to visualize and are plotted as a series of connected vertical and horizontal dashed lines between the Constraint Method tradeoff points. Obviously, they do not provide any insightful information. In contrast, even in the absence of assuming convexity, valuable auxiliary information is available for the bounds of the tradeoff curve if the TSR Algorithm is used. Clearly, the TSR Algorithm provides much more information as to the location of the tradeoff curve between the tradeoff points when compared to the Constraint Method. The TSR partial tradeoff curve boundaries provide information to the decision-maker or analyst that can be used to determine if more tradeoff points should be generated in between the present tradeoff points. For example, if the decision-maker is contemplating spending approximately 50 instead of 40 million $US/yr in order to improve the water quality system reliability, the decision-maker would know from Figure 5.11 that at a cost of 50 million $US/yr the maximum possible reliability will be approximately 0.85. If this maximum possible reliability value doesriotsatisfy the decision-maker's minimum desired reliability level for the proposed spending increase, the decision-maker would be able to reject the increase in spending without evaluating another 116 tradeoff point. This is a significant advantage of the TSR Algorithm, especially for computationally intensive problems like the one considered in this work. Figure 5.11. Cost-reliability tradeoff curves and partial tradeoff boundaries for the TSR and Constraint Methods Figure 5.12 shows the same tradeoff curve that appears in Figure 5.11 for the TSR Algorithm and a new tradeoff curve for the Constraint Method that is developed with the constrained reliability levels specified to correspond to the reliability values of the optimal tradeoff points from the TSR Algorithm. The number of GA generations required to find each interior tradeoff point is also shown for each tradeoff curve. Although both methods produce non-inferior tradeoff points and essentially the same tradeoff curves, the TSR Algorithm produces the tradeoff curve more efficiently in terms of the number of required GA generations. An average of 287 GA generations are required to find the optimal tradeoff points using the TSR Algorithm as compared to an average of 559 GA generations required to find the optimal tradeoff points using the Constraint Method. This efficiency difference is due to the use of a penalty function with the GA. The Constraint Method relies heavily on the penalty function to guide the search to 117 solutions with acceptable reliability levels and therefore reduces the efficiency of the GA relative to the nearly unconstrained search when the TSR Algorithm is used. The TSR Algorithm, as applied here, only utilizes a penalty method to ensure solutions between the reference points are considered and thus the penalty function values are rarely greater than zero. Figure 5.12. TSR and Constraint Method comparison at approximately the same reliability levels The TSR Algorithm is shown to produce a more representative tradeoff curve than the Constraint Method if the same number of tradeoff points is generated by each technique. In addition, in the GA implementation of the TSR Algorithm and Constraint Method, the GA converges faster under the TSR Algorithm than under the Constraint Method. Therefore, the TSR Algorithm is used as the multi-objective method to generate the remaining tradeoff curves in this case study. 118 5.2.3 FORM and MCS Comparison in the Cost-Reliability TSR Model FORM and MCS are compared as the reliability estimation technique in the Cost-Reliability TSR model (Equation (34)). All FORM and MCS tradeoff curves in this section are generated such that the total computational time required for all tradeoff points is approximately the same. The minimum number of FORM performance function evaluations is calculated as 39 from Equation (25) when six random variables are used to characterize the system uncertainty. Therefore, two tradeoff curves using 40 and 80 MCS realizations to estimate the reliability of each solution generated by the GA are developed for comparison with the FORM-based TSR cost-reliability tradeoff curve. When 40 and 80 MCS realizations are used as the reliability estimation technique within the water quality management model, 1000 and 500 GA generations respectively, are used to generate each tradeoff point. The FORM and MCS reliabilities estimated within the water quality management models are referred to as the nominal reliability estimates while the more accurate reliability, estimated using a consistent set of 5000 MCS realizations, for the optimal solutions is referred to as the actual reliability estimate for that solution. The actual reliability estimates are presented for each cost-reliability tradeoff curve. For reasons described in Section 4.3.3.1 all MCS tradeoff points require post-processing the GA output to find the optimal solution. For each MCS tradeoff point this post-processing includes identifying the best two or three solutions that have a reliability estimate accumulated by solutions that persist as the elite population member over a sufficient number of generations and then determining the actual reliability estimate of these two or three solutions. A sufficient number of generations for each tradeoff point is roughly determined by the number of MCS realizations required to achieve approximately the same reliability accuracy (using Equation (26) at the 95% confidence level) as the FORM-based reliability estimates in the same region of the tradeoff curve. Comparing the FORM and actual reliability estimate determines the FORM reliability error of each optimal solution on the tradeoff curve. The solution with the best fitness when considering the actual reliability is then deemed the optimal MCS solution. Figure 5.13 compares the FORM- and 40-realization MCS-based TSR cost-reliability tradeoff curves with the nominal and actual reliability estimates. The solid lines represent the actual reliability tradeoff curves while the dashed lines represent the nominal tradeoff curves. The reference points for the FORM-based nominal TSR tradeoff curve are defined by the FORM 119 reliability estimates while the reference points for the MCS-based nominal TSR tradeoff curve are based on the MCS reliabilities estimated with 5000 realizations. This means that MCS and FORM are being utilized for slightly different optimization problems. The FORM nominal estimates of reliability underestimate the actual reliability by an average of 0.052. This is consistent with the results in Section 5.1.1 except that the FORM reliability error is now slightly higher. The size of the FORM reliability error suggests that using a small number of MCS realizations to roughly estimate each solution's reliability may also be a comparable method for reliability estimation within an optimization model. Comparing the nominal and actual MCS reliabilities shows that the error in the nominal MCS reliability estimate is equal to an average value of 0.017, which is significantly less than the FORM reliability error. In addition, the nominal MCS tradeoff curve suggests initially that one of the nominal FORM tradeoff points is grossly inferior. However, comparing the actual FORM and MCS tradeoff curves shows that neither reliability estimation technique produces inferior tradeoff points. The differences between the nominal and actual cost-reliability tradeoff curves demonstrate the need to do a postoptimization check on all optimal solution nominal reliability estimates with a more accurate MCS estimate of reliability. 120 20 -I 0.3 , , , , , , | 0.4 0.5 0.6 0.7 0.8 0.9 1 Reliability Figure 5.13. TSR Algorithm tradeoff curves using FORM and 40 MCS realizations for reliability estimation 120 Focusing on only the FORM and MCS actual tradeoff curves in Figure 5.13, it is clear that since the FORM actual tradeoff curve points are furthest from the reference lines used in either the FORM or MCS nominal tradeoff curves, the use of FORM as the reliability estimation technique in the water quality management model produces a better representation of the tradeoff curve. In terms of computational efficiency, the optimal FORM solutions are found in an average of 287 GA generations while the optimal 40-realization MCS solutions are found in an average of 652 generations. Therefore, the FORM-based tradeoff points are found over twice as quickly as the MCS-based tradeoff points. In terms of actual computer processor time, 1000 GA generations for each FORM and MCS tradeoff point requires about 19.5 and 20 hours respectively, on a Pentium III 600 MHz personal computer. Figure 5.14 shows the FORM, 40-realization MCS and 80-realization MCS tradeoff curves using the actual reliability estimates. Since the MCS and FORM extreme tradeoff curve points are the same as in Figure 5.13, only the three interior points in each tradeoff curve are shown. The 80realization MCS tradeoff curve points are much closer to the FORM tradeoff curve points than the 40-realization MCS tradeoff curve. However, the 80-realization MCS tradeoff point with the highest cost and reliability is inferior to the corresponding FORM tradeoff point. An average of 313 GA generations is required to converge to the 80-realization MCS tradeoff points. Since 80 MCS realizations is approximately double the minimum number of performance function evaluations required by FORM and FORM is observed to most often converge with a minimum number of performance function evaluations, 313 GA generations of an 80-realization MCS tradeoff point requires approximately the same computational time as 626 GA generations of a FORM-based tradeoff point. Therefore, compared to the FORM tradeoff points, which took an average of 287 GA generations to converge, more than twice as much computational time is required for the 80-realization MCS tradeoff points. The comparison of FORM and MCS in the TSR water quality management model suggests that when FORM is used as the reliability estimation technique, it produces the most representative and efficient tradeoff curve. This result is most likely due to the observed consistency of the FORM reliability error as compared to the random 40- or 80-realization MCS reliability error. All FORM and 5000-realization MCS reliability comparisons in this research showed that FORM always over-predicts the actual reliability. This means that even though the accuracy of the FORM-based reliabilities may be somewhat poor, the solution that is truly the most reliable 121 is likely to be identified as such by the FORM-based reliability estimates. Similarly, the difference between actual reliability values for any two solutions is very likely to be closely replicated in the nominal FORM reliabilities for these solutions. In comparison, the 40- or 80realization MCS reliability error relative to the actual reliability estimates will not be consistent for various solutions due to the random generation of each set of MCS realizations per GA generation. For example, in one generation, a set of MCS realizations could determine that solution A is more reliable than solution B while a new set of MCS realizations in the next generation could determine that solution B is more reliable than solution A. This randomness prevents the MCS-based water quality management models from consistently improving the fitness throughout the GA generations and instead produces, for the elite solutions, a fitness through generations that resembles random noise with no significant trend. As a result, good solutions are not being selected for reproduction throughout all generations and valuable genetic information is continually lost from the population. 122 Despite the above phenomenon, Figure 5.14 shows that the MCS-based tradeoff curve with 80 realizations is reasonably close to the FORM-based tradeoff curve. Further improvements in the use of MCS within a GA could greatly reduce the loss of genetic information due to the random MCS reliability errors and possibly produce an even closer or better cost-reliability tradeoff curve than the FORM-based tradeoff curve. Changes in the GA that would accomplish this include the accumulation of the nominal reliabilities during the optimization and the storage of elite solutions from earlier generations for future reinsertion to the population if the later generation elite solutions were eventually identified as having smaller fitness values than the elite solutions they replaced. Another improvement would be the determination of an optimal number of MCS realizations for nominal reliability estimation within the optimization model. Since these GA-MCS improvements are not within the scope of this research, FORM is used to estimate the nominal values of the performance indicators within the water quality management models for all further tradeoff curves. 5.2.4 Tradeoff Curve Comparison Between Cost and Performance Indicators FORM-based cost-vulnerability and cost-resilience tradeoff curves are generated using the TSR Algorithm. Due to the increased FORM computational burden of vulnerability and resilience estimation relative to reliability estimation, the maximum number of GA generations in the water quality management models is reduced from 1000 to 350. Again, the optimal solutions for all performance indicators are checked with 5000 MCS realizations to achieve a more accurate estimate of the actual performance indicator values. Tradeoff curves plotted with these MCS estimates of the performance indicators are referred to as the MCS checked cost-vulnerability or cost-resilience curves. Since the TSR partial tradeoff curve bounds only apply to the nominal FORM estimates of vulnerability and resilience instead of the more accurate MCS checked performance indicator values, no tradeoff curve bounds are plotted for the cost-resilience and cost-vulnerability tradeoff curves. In other words, bounds on the tradeoff curve can be constructed when the FORM vulnerability and resilience values are used, but they cannot be constructed with the MCS estimates of the optimal solutions vulnerability or resilience value since the TSR Algorithm did not use the MCS estimates. All optimal solutions have their two, yet to be estimated, performance indicator values assessed outside of the water quality management models by FORM and 5000 MCS realizations. 123 Figure 5.15 shows the FORM and MCS checked cost-vulnerability tradeoffs as generated by the TSR Algorithm. The FORM vulnerability errors are similar to those observed in Section 5.1.1. The small differences in the information provided by the cost-reliability and cost-vulnerability tradeoff curves are due to the additional consideration of the water quality system failing with a DO level less than 3 mg/L. Since this failure probability is quite small, the information given by the cost-reliability and cost-vulnerability curves is essentially the same. The average computer processor time that is required for a single cost-vulnerability tradeoff point is approximately 13.5 hours on a Pentium III 600 MHz personal computer. Consistent with the fact that the vulnerability estimation requires FORM to estimate two reliability values, the processor time that is required per GA generation for the cost-vulnerability tradeoff points is approximately double that observed for producing points on the cost-reliability tradeoff curve. 120 0.5 1.0 1.5 2.0 2.5 Vulnerability Figure 5.15. TSR cost-vulnerability tradeoff curves for the FORM and 5000-realization MCS estimates of vulnerability 124 Figure 5.16 shows the FORM and MCS checked cost-resilience TSR tradeoff curves. In addition, the approximate 95% confidence limits for the actual resilience value are also shown for the three interior tradeoff points. The MCS checked tradeoff curve is very flat from the lowest resilience value until a resilience of approximately 0.87. This means that the water quality system resilience can be improved for small increases in total cost up to a resilience of approximately 0.87. Any further improvement in resilience above 0.87 requires a large increase in total waste treatment costs. A i â€” A â€” F O R M tradeoff - - -A - â€¢ M C S checked tradeoff - - - Approximate 95% confidence limits on M C S resilience I & < k 0.4 0.5 * t A 0.6 0.7 Resilience 0.8 0.9 1.0 Figure 5.16. T S R cost-resilience tradeoff curves for the F O R M and 5000-realization M C S estimates of resilience Referring to the upper portion of the MCS checked cost-resilience tradeoff curve in Figure 5.16, it appears that the water quality management model solution at a cost of 60 million $US/yr may be inferior. If this is in fact true then this result may indicate that FORM is not suitable for resilience estimation in this case study. However, out of the 5000 MCS realizations utilized to check each of the optimal solution resilience estimates, only 86 to 225 resulted in failure in the first time step. Therefore, the accuracy of the MCS estimated resilience values is based on 86 to 125 225 realizations rather than 5000 realizations and therefore greatly reduces the quality of the MCS resilience estimates for the optimal solutions. The approximate 95% confidence limits on the actual resilience value show that it cannot be concluded that the FORM-based TSR Algorithm has identified an inferior tradeoff point. The inclusion of two of the three FORM resilience estimates within the approximate 95% confidence limits for the actual resilience value is an improvement over the FORM resilience accuracy displayed in Figure 5.3(a). The computational time required for each GA generation in each cost-resilience tradeoff point must be a minimum of approximately four times more than the computational time required for each GA generation in each cost-reliability tradeoff point based on Equations (25) and (27). This is observed to be the case for one of the resilience points where 350 generations required approximately 34 hours of processor time on a Pentium III 600 MHz personal computer. However, for the other two tradeoff points, FORM fails to converge to a solution in at least one resilience evaluation and thus pauses the execution of the water quality management model program. For one of these points, FORM fails to converge during a number of resilience calculations and results in the optimization trial requiring more than 5 days to complete. Not all of this time can be attributed to processor time since the program often pauses. However, it is still observed that the actual processor time is still greater than 34 hours. When FORM does not converge, the solution is assigned a poor resilience estimate. The impact of this on the optimization results is not quantitatively assessed but this arbitrarily assigned resilience value could potentially result in a good solution being lost from the population due to FORM nonconvergence. The inability of FORM to converge for a number of candidate solutions in the cost-resilience water quality model is a definite disadvantage of the FORM technique when used to estimate the performance indicators in an optimization model. 126 Table 5.3 summarizes the FORM and 5000-realization MCS estimates of the performance indicators for the cost-reliability, -vulnerability and -resilience tradeoff curve solutions. For example, when reliability is the performance indicator in the water quality management model and tradeoff point 3 is considered, the vulnerability and resilience estimates by FORM and MCS correspond to the cost-reliability optimal solution at tradeoff point 3. Comparing solutions with similar total costs shows that none of the identified solutions are inferior to each other and that the performance indicators may not be directly conflicting. In other words, if costs are slightly reduced, then reliability and resilience levels decrease and vulnerability increases. However, it cannot be concluded from the results in Table 5.3 that an optimal cost-reliability solution will also have an optimal resilience or vulnerability for that particular cost, or vice-versa. In order to determine this, the three performance indicators would have to be optimized individually with costs constrained at the same level. The computational time required to assess the performance indicator values of the optimal solutions using 5000 MCS realizations is small in comparison to the optimization trials. MCS reliability and vulnerability assessment requires approximately 0.5 hours while resilience requires between 0.5 and 0.75 hours, depending on the failure probability in the first time step, on a Pentium III 600 MHz personal computer. 127 c o o fi o U 00 r-; 00 CN o "cu H ON cn cn CN ON 00 CO CO ON o CO ON CO CO CO NO NO o NO NO NO NO CO CO CO ON r- CO Â© CO CO o Â© CO NO CD s CU X! -O CD 00 c U _o -4â€”* u o â€¢ EH o â€¢3 a Pi rfi o co o es u ccj _N CD cu u O es NO o NO ON w a IT) IT) oo o NO ON ON 00 oo oo r H ON NO 00 ON NO O oo CN 00 oo NO OO 00 CO ON oo NO O r H o o S .2 NO ON oo ON o CO NO o ON CN NO ON NO NO Â°0 0O Â°0 co ON C5 OO oo Â£ -5 t. o Q. a> o 5 s oo oo < ~ Â£ => co <U O ON CO NO -3NO NO CN CN CN CN O Â© CN cn ON NO ON o O o o CN IT) NO C\ VO CO co NO CU CN NO CN CO CN Â© 'es' T3 cs NO o 00 CN o 5 o o â€¢. o N CD O oo IT) CO NO OO CO ON OO IT) IO O Â© CN >>U NO CN o CN CO oo o NO NO -J- r H ON ON 00 NO NO O O 0O oo oo NO co CO o oo 00 o CN NO NO ON in IcO CN o ON CN IT) OO Â£ Â£ CO ON Â© Â© s !eS ON ^H CN CO H CN CO CN CO fi a 3 3 C Z J CD o fi CU cs 1-1 CD H PH 1 > CD co CD Pi 'rH CN 128 6 CONCLUSIONS AND FUTURE WORK This research demonstrates a number of novel approaches and techniques in the multi-objective management of water quality systems under uncertainty. Some of the conclusions of this work apply specifically to the Willamette River case study while others apply to the general water quality management framework developed. Many of the findings are applicable across all disciplines that attempt to model systems with uncertainty or manage systems with multiple objectives. In addition, this work highlights areas where future research is needed. FORM and MCS is applied to the Willamette River Basin case study in order to estimate the water quality system reliability, vulnerability and resilience with respect to meeting various DO standards at the end of the river. To the author's best knowledge, this, and the author's work in Maier et al. (1999), is the first application in water resources where FORM has been used for vulnerability and resilience estimation. The results clearly show that most of the FORM estimates of the performance indicators are significantly different than the corresponding MCS estimates. However, when the relative efficiencies of FORM and MCS are assessed assuming both methods are required to achieve the same accuracy level for the performance indicator estimates, FORM appears to be the most efficient performance indicator estimation technique. Also, auto-correlation of the water quality model Ka and SOD uncertainty sources are shown to greatly impact the system resilience. This exemplifies the need to have reasonable instead of arbitrary estimates of the auto-correlations for the Ka and SOD uncertainty sources. The FORM and MCS efficiency comparison method used here is unique from most FORM and MCS reliability comparisons in the literature. FORM is most often compared with MCS probabilities that are based on an arbitrarily large number of MCS realizations or a sufficient number of realizations that reach some level of convergence in the probability estimate (Tung 1990; Melching 1992; Melching and Anmangandla 1992; Maier et al. 1999). The FORM and MCS efficiency studies by Melching and Anmangandla (1992) and Maier et al. (1999) conclude that FORM reasonably approximates MCS but that FORM is much more efficient than MCS. This conclusion is made without considering that MCS will most often require fewer realizations to achieve the observed FORM accuracy level than the realizations required for MCS to converge. 129 A GA is utilized as the optimization technique to generate points on two-objective tradeoff curves with combinations of total annual CBOD waste treatment costs and the water quality system reliability, vulnerability or resilience as the objectives for the Willamette River case study. To the author's best knowledge, this, and the author's work in Vasquez et al. (2000) are the only two examples of GAs used in conjunction with FORM to estimate tradeoffs involving waste treatment costs and any of the performance indicators in water resources. The practically unrestricted applicability of GAs is highlighted in this work by the ability of the GA to operate freely from the complex determination of the performance indicator estimates for each trial solution. In addition, the ability of the GA to use the actual discrete waste treatment options as decision variables is very desirable since this does not require any additional approximation or linearization of the waste treatment cost data. The disadvantage of the GA optimization approach is also highlighted by the large number of fitness function evaluations required and the resultant computational burden that is on the order of days for tradeoff curve generation. Although this is a large computational burden, it would seem unreasonable to reject such an approach in the real design of multi-million dollar wasteload allocation programs solely on the basis that the computational time would be a few days. Further, parallel computation can be used to greatly improve GA computational time. The FORM and MCS comparison is also extended to their respective inclusion as the performance indicator estimation technique in the water quality management models. FORM produces better tradeoff curves than MCS with either 40 or 80 MCS realizations per reliability estimate. The FORM tradeoff points are also found in a fewer number of GA generations than the MCS tradeoff points. However, the 80-realization MCS tradeoff curve resembles the FORM tradeoff curve closely. Besides FORM outperforming MCS in terms of the tradeoff curves developed, the FORM and MCS comparisons identified two additional positive aspects regarding FORM for the estimation of the performance indicators. FORM efficiency relative to MCS is highest at higher reliability and resilience values - on the order of 10 to 50 times as efficient as MCS. This is a significant advantage of FORM because most water quality systems are managed to achieve high values of reliability or resilience. Another observed benefit of FORM is the consistency of the FORM performance indicator error, especially in reliability estimation. The consistent FORM reliability error judges the relative reliability of various solutions correctly and is therefore the primary 130 reason that FORM-based cost-reliability tradeoffs are better than the MCS-based tradeoffs. While the former benefit is quite possibly a result that can be expected in other case studies, the later should not be taken as a general result. For example, if in another case study, the management decisions included locating a number of point sources along a river and/or flow regulation levels, in addition to waste treatment levels at point sources, the likelihood of FORM exhibiting a consistent reliability error over all management decisions is small. The list of disadvantages of the FORM technique is more extensive than its positive aspects. Some of these disadvantages include the decrease in FORM efficiency with an increase in random variables, non-convergence of FORM in some cases, the possible identification of drastically different design points for very similar solutions and the fact that FORM accuracy is problem specific. In contrast, when MCS is used to estimate the performance indicators, the efficiency is independent of the number of random variables, MCS will always give a point estimate of the performance indicator, MCS samples over the entire random variable space and thus multiple design points are accounted for in all MCS estimates and MCS accuracy is only dependent on the number of realizations and the actual probability value it is approximating. The most significant disadvantage of FORM listed above is the efficiency dependence on the number of random variables. When developing an uncertain model for any system, one would ideally like to be unconstrained in the determination of the sources of uncertainty in that the uncertainty estimation method should not dictate or greatly influence the number or type of random variables selected. This is especially true considering that the work by Thomann (1998) suggests the future of water quality modelling lies in massive ecosystem based models where there will be an abundance, perhaps even hundreds, of significant sources of uncertainty. The other big disadvantage of the FORM technique is that it, as applied here, is not applicable to a system that is modelled dynamically. Resilience is estimated in this research using a steadystate model (see Section 4.2.2). Since most water quality or water resource systems where resilience is an important management objective will likely be modelled with dynamically, the FORM resilience approach presented here is not applicable. In contrast, the MCS resilience approach could still be valid for resilience estimation using a dynamic model. For example, one MCS realization would be the generation of a single set of input time series for the uncertain variables, rather than only two time steps of inputs. 131 Despite the disadvantages of FORM, it is still shown to produce reasonable cost-vulnerability and cost-resilience tradeoff curves as the performance indicator estimation technique in the water quality management models. The tradeoff curves are reasonable in that none of the optimal solutions are inferior to one another. Since MCS based comparisons with the optimal FORM performance indicator estimates require a relatively small amount of computation time and significant inaccuracies are observed in the optimal FORM performance indicator estimates, checking all FORM estimates is an important post-optimization step in the water quality management process. This is especially true if the inaccuracies of the FORM performance indicator estimates cover a correspondingly large range of management costs. Some previous studies (e.g. Xu and Goulter 1999) that determine a cost-reliability tradeoff with an approximate reliability estimation technique have failed to check any of the optimal reliability estimates with a large number of MCS realizations. Overall, FORM seems to be a valid approach for estimation of reliability, vulnerability and resilience in the Willamette River Basin case study. However, it cannot be concluded that FORM is a better technique for performance indicator estimation than MCS within an optimization model. To more clearly determine the relative merits of the FORM and MCS performance indicator estimation techniques the following topics should be investigated further in future work: â€¢ The FORM and MCS efficiency comparison here should be reassessed in terms of the confidence level used to estimate the number of required MCS realizations since the comparisons here require that approximately 95% of MCS estimates are at least as or more accurate than the FORM-based values of reliability and resilience. For example, it could be argued that MCS is a more efficient method than FORM if, for the same number of performance function evaluations, MCS performed more accurately than FORM in just over 50% of MCS trials. â€¢ Other techniques that significantly improve MCS or FORM efficiency or accuracy, such as Importance Sampling, should be used to estimate reliability, vulnerability and resilience in further FORM and MCS comparisons. â€¢ Improvements to the combination of GAs with MCS for optimization under uncertainty should be undertaken in order to fairly compare the optimization results that use FORM and MCS as the performance indicator estimation technique. 132 â€¢ A two-step GA optimization method should be undertaken in which FORM is used to quickly identify good solutions and then the GA is restarted from the good FORM solutions using MCS to refine the solution with more accurate estimates of the performance indicators. â€¢ Alternative measures of system vulnerability should be defined in order to characterize the system vulnerability by a measure that has physical meaning and better test FORM for vulnerability estimation. For example, if the failure regions are defined in terms of the spatial extent of the DO concentration being less than the standard and the corresponding vulnerability weights represent the spatial extent of the failure events, then the vulnerability measure can be interpreted as the expected length of the river failing to meet the DO standard. FORM'S ability to estimate the required failure probabilities would then need to be tested since the performance function in this case would include the length of the river in failure instead of a point estimate of the DO concentration. â€¢ The multi-objective water quality management models developed for the Willamette River should be expanded to include the simultaneous management of all significant pollutants such as nitrogen and phosphorous. Recent work by Carmichael and Strzepek (2000) demonstrates the need for inclusion of interactive pollutant effects in multiple-pollutant deterministic wastewater management by showing that for their case study managing only BOD emissions leads to non-optimal solutions when nitrogen and phosphorous emissions are considered. It is reasonable to expect this result is also true for stochastic wastewater management. This research has also presented the TSR Algorithm as a potentially novel approach to approximating convex tradeoff curves in the most efficient and representative way possible for multi-objective problems. The TSR Algorithm can be applied to multi-objective problems in any field with any number of objectives, used as an interactive or non-interactive method with a decision-maker and implemented using any optimization technique that is valid for the objectives considered. The TSR Algorithm, while similar to other multi-objective approaches such as Das (1999) and Fandel (1972), is novel because it is used to determine the entire tradeoff curve, provides a maximum amount of auxiliary information on the bounds of the tradeoff curve and produces a mean estimate of the tradeoff curve. Furthermore, the mathematical formulation of the TSR Algorithm contains one less decision variable than the approach proposed by Das (1999) forfindingthe knee of the tradeoff curve. 133 Previous work has highlighted the problems of the Weighting Method (e.g Das and Dennis 1997) and the results in this work show that the TSR Algorithm produces a much more representative tradeoff curve than the Constraint Method when each method generates the same number of tradeoff points. Therefore, the TSR Algorithm is a superior algorithm for generating tradeoff curves in convex two-objective problems when compared to the commonly used the Weighting or Constraint Methods. The case study also shows that the TSR Algorithm provides valuable information regarding the boundaries of the tradeoff curve location even if convexity of the tradeoff surface is not assumed. The method can also identify non-convex or partially nonconvex regions of the tradeoff curve where decision-makers or analysts can then decide if additional tradeoff points need to be ascertained by another multi-objective method. Future work regarding the TSR Algorithm should include the following: â€¢ Application of the TSR Algorithm to multi-objective problems with more than two objectives to determine if any minor modifications to the algorithm are necessary for application in with more than two objectives. This should be followed by further comparisons with the NISE Method (Cohon et al. 1979) and its derivative method in Solanki et al. (1993) and then a further review of the multi-objective literature to determine if the TSR Algorithm is in fact a novel multi-objective approach. â€¢ Demonstration of the TSR Algorithm using other optimization methods such as linear and nonlinear programming for well known multi-objective examples in the literature. â€¢ Comparison of the TSR Algorithm with MOGAs to compare tradeoff surfaces generated and relative efficiencies. â€¢ Combination of the TSR Algorithm and MOGAs to map out tradeoff surfaces in the most efficient manner possible. MOGAs strive to generate an even spread of points on the tradeoff surface, even if there is little tradeoff observed between objectives. This behaviour finds relatively useless solutions and adds unnecessary computational time since in most management situations, decisions will not be made in such level areas of the tradeoff surface because this implies that one or more objectives can be improved with little cost to the other objectives. Therefore, a logical and efficient approach to multi-objective problems, especially for case studies like this that are very computationally intensive, would be to generate the tradeoff surface in two steps. First, the TSR Algorithm would be used to get a quick and accurate representation of the tradeoff surface. Then, for areas of the tradeoff 134 surface where the decision-maker decides he or she is most interested a MOGA could be used to find a variety of evenly spread candidate solutions. Similarly, this could also be the approach for partially non-convex problems where the TSR Algorithm first identifies the non-convex regions and then, if the non-convex regions are of interest, a MOGA could be used to map out the tradeoff curve in the respective regions. Application of the TSR Algorithm, or related variants, outside of the multi-objective optimization field for the approximation of any continuous functions as piecewise linear curves. For example, in the FORM approach to reliability estimation, the failure surface is approximated in standard normal space by a hyperplane tangent to the design point (see Figure 3.1) and the objective of FORM is to approximate the area of the failure domain that is bounded the failure surface. The TSR algorithm can be iteratively applied in standard normal random variable space to maximize the perpendicular distance between a hyperplane, defined by some predefined reference points that are identified initially by the minimum and maximum values of the random variables, and the failure surface. With a suitable optimization technique constrained to only search the failure surface between the reference points, the TSR Algorithm or a closely related variant may be able to quickly and efficiently find a piecewise linear surface that closely approximates the failure surface. This modified FORM technique would generate more accurate reliability estimates than the traditional FORM method and would even provide an estimate of the maximum error in each reliability estimate based on the construction of the boundaries on the failure surface location. Furthermore, the TSR Algorithm does not require a convexity assumption to be applied in standard normal random variable space since the failure surface can be found on either side of the reference hyperplane. Therefore, the TSR Algorithm may be even better suited for approximation of the failure surface in reliability analysis than multi-objective tradeoff surface approximation. 135 NOTATION The following symbols are used in this thesis: A = penalty multiplier; B = penalty exponent; c = a constant; D () = CDF of the load, L; dj = normalized scalar distance in objective i from the reference hyperplane to the L value of the ith objective; E[ ] = expected value; ej = probability that the system performance resides in discrete failure state j; F = failure domain; fx(x) = joint probability density function (PDF) of A!"; G(X) = performance function, where the failure (limit state) surface G = 0 separates all combinations of X that lie in the failure domain from those in the survival domain; g = perpendicular distance in objective space, in the direction of the ideal point, from the vector of objective function values (Z(y)) to the reference hyperplane; ' h (Y,Q ,X ) = estimated reliability for a given vector of discrete waste treatment options, s par water quality standard at the critical location in the stream and vector of probability distribution parameters describing the random variables; i = general index; j = general index; k (Y) = total cost of waste treatment; L = a system's load m = number of objectives in a multi-objective problem; NF = number of performance function evaluations required by FORM for a reliability estimate; NFres = number of performance function evaluations required by FORM for a resilience estimate; NMCS = number of MCS realizations used to estimate reliability; 136 NMCS cond - number of times the second performance function (for time step two) is conditionally evaluated for a MCS resilience estimate; NMCS res - total number of performance function evaluations required for a MCS resilience estimate; NMcsresi total number of MCS realizations used in the first time step for a MCS = resilience estimate; n â€” number of random variables in X for a reliability or vulnerability estimate; n â€” number of random variables in X for a resilience estimate; Pf\ = probability of failure for failure mode 1; pj\2 = joint probability of failure for failure modes 1 and 2; Q = water quality standard at a critical location; q = number of FORM iterations; R = a system's resistance; r = (r\, n, â€¢â€¢â€¢ n) is a vector of constant coefficients in the equation of a hyperplane; r = slope of a line; S = survival domain; s = example value of a decision variable; Tf = length of time a system's output remains unsatisfactory after a failure; t = index for a time step; V = vector of uncorrected standard normal variables transformed from the res s correlated random variables that influences a system's load and resistance; V* = realization of the vector of uncorrelated standard normal variables at the design point; Wj = numerical indicator of the severity of failure state j; X = vector of random variables that influences a system's load and resistance; X* = realization of the vector of random variables at the design point, or the most likely failure point, the point on the failure surface closest to the mean point; X = vector of probability distribution parameters describing the random variables; Y - Y = decision variable space; y = vector of decision variable values; Zj = management objective i; par s Zj ax m = vector of discrete waste treatment options; upper bound of minimization objective i; 137 Z*i = minimum value for objective i; Z* = ideal point in objective space for two objectives that are bounded by points A AB and B; Z = objective space; Z(y) = vector function of objectives; Z = (Zi, Z 2 , Z â€ž ) s T is a vector of the unknown variables in the hyperplane equation and is also a vector of the objectives; = standard normal variate with cumulative probability level of OJ/2; a = reliability; a = B = reliability index; cp( ) = standard normal CDF; (p( , ;p) = CDF for a bivariate normal vector with zero mean values, unit variances and s specified reliability level for the system; correlation coefficient p\ y - v = vulnerability; cos 61 = direction cosine of random variable 1 in the standard normal space; 0 = angle between the vector of the design point and the axis of the random n resilience; variable Xâ€ž, and <74 , ;p) â€” PDF for a bivariate normal vector with zero mean values, unit variances and correlation coefficient p. 138 BIBLIOGRAPHY Aly, A. H., and Peralta, R. C. (1999a). "Comparison of a genetic algorithm and mathematical programming to the design of groundwater cleanup systems." Water Resources Research, 35(8), 2415-2425. Aly, A. H., and Peralta, R. C. (1999b). "Optimal design of aquifer cleanup systems under uncertainty using a neural network and a genetic algorithm." Water Resources Research, 35(8), 2523-2532. Beck, M. B. (1987). "Water Quality Modeling: A Review of the Analysis of Uncertainty." Water Resource Research, 23(8), 1393-1442. Brown, L. C , and Barnwell, T. O., Jr. (1987). "The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: documentation and user manual." Report EPA/600/387/007, US Environmental Protection Agency, Athens, GA. Burges, S. J., and Lettenmaier, D. P. (1975). "Probabilistic methods in stream quality management." Water Resources Bulletin, 11(1), 115-130. Burn, D. H., Venema, H. D., and Simonovic, S. P. (1991). "Risk-based performance criteria for real-time reservoir operation." Canadian Journal of Civil Engineering, 18(1), 36-42. Burn, D. H., and Yulianti, J. (1999). "Genetic algorithms for managing water quality systems under uncertainty." Journal of Water Resources Planning and Management, submitted. . Cardwell, H., and Ellis, H. (1993). "Stochastic dynamic programming models for water quality management." Water Resources Research, 29(4), 803-813. Carmichael, J. J., and Strzepek, K. M. (2000). "A multiple-organic-pollutant simulation/optimization model of industrial and municipal wastewater loading to a riverine environment." Water Resources Research, 36(5), 1325-1332. Chadderton, R. A., Miller, A. C , and McDonnell, A. J. (1982). "Uncertainty analysis of dissolved oxygen model." Journal of the Environmental Engineering Division, 108(EE5), 1003-1013. Chan-Yan, D. (2000). "Reservoir turbidity: Neural network modelling and the evaluation of performance indicators," M.A.Sc, University of British Columbia, Vancouver. Chen, H. W., and Chang, N.-B. (1998). "Water pollution control in the river basin by fuzzy genetic algorithm-based multiobjective programming modeling." Water Science & Technology, 37(3), 55-63. Chen, Y. M. (1997). "Management of water resources using improved genetic algorithms." Computers and Electronics in Agriculture, 18(2-3), 117-127. Cheng, S.-T. (1982). "Overtopping risk evaluation for an existing dam," Ph.D., University of Illinois at Urbana-Champaign. Cieniawski, S. E., Eheart, J. W., and Ranjithan, S. (1995). "Using genetic algorithms to solve a multiobjective groundwater monitoring problem." Water Resources Research, 31(2), 399-409. Cohon, J. (1978). Multiobjective programming and planning, Academic Press, New York. Cohon, J. L., Church, R. L., and Sheer, D. P. (1979). "Generating multiobjective trade-offs: an algorithm for bicriterion problems." Water Resources Research, 15(5), 1001-1010. Correia, F. N., Santos, M. A., and Rodrigues, R. R. "Risk, resilience and vulnerability in regional drought studies." Systems Analysis Applied to Water and Related Land Resources, Proceedings of the IFAC Conference, Lisbon, Portugal, 139-144. Dandy, G. C , Simpson, A. R., and Murphy, L. J. (1996). "An improved genetic algorithm for pipe network optimization." Water Resources Research, 32(2), 449-458. Das, I. (1999). "On characterizing the "knee" of the Pareto curve based on normal-boundary intersection." Structural Optimization, 18, 107-115. Das, I., and Dennis, J. (1998). "Normal-boundary intersection: a new method for generating the Pareto surface in nonlinear multicriteria optimization problems." Society for Industrial & Applied Mathematics Journal of Optimization, 8(3), 631-657. Das, I., and Dennis, J. E. (1997). "A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems." Structural Optimization, 14, 63-69. Der Kiureghiam, A., and Liu, P.-L. (1986). "Structural reliability under incomplete probability information." Journal of Engineering Mechanics of the American Society of Civil Engineers, 112(1), 85-104. Duckstein, L., and Bernier, J. "System framework for engineering risk analysis." Risk-Based Decision Making in Water Resources, Proceedings of an Engineering Foundation Conference, Santa Barbara, CA, USA, 90-110. Edward C. Jordan Co., I. (1979). "Preliminary data base for review of BATEA effluent limitations guidelines, NSPS, & pretreatment standards for the pulp, paper, & paperboard point source category." NOTE: not an official Environmental Protection Agency publication, U.S. Environmental Protection Agency, Portland, Maine. 140 Eschenauer, H., Koski, J., Osyczka, A., and (Editors). (1990). Multicriteria design optimization, Springer-Verlag, New York. Fandel, G. (1972). Optimale Entscheidung bei mehrfacher Zielsetzung, Springer Verlag, Berlin. Fiering, M. B. (1969). "Forecasts with varying reliability." Journal of the Sanitary Engineering Division of the American Society of Civil Engineers, 95(4), 629-644. Fiering, M. B. (1976). "The role of systems analysis in water program development." Natural Resources Journal, 16. Foschi, R. O., and Folz, B. (1990). "RELAN: RELiability ANalysis." , University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z4. Fuessle, R. W. (1987). "Air quality planning: A general chance-constraint model." Journal of Environmental Engineering, 113(1), 106-123. Goldberg, D. E. (1989). Genetic algorithms in search, optimization and machine learning, Addison-Wesley, Reading, MA. Goldberg, D. E., and Kuo, C. H. (1987). "Genetic algorithms in pipeline optimization." Journal of Computing in Civil Engineering, 1(2), 128-141. Halhal, D., Walters, G. A., Ouazar, D., and Savic, D. A. (1997). "Water network rehabilitation with structured messy genetic algorithms." Journal of Water Resources Planning and Management, 123(3), 137-146. Hashimoto, T., Stedinger, J. R., and Loucks, D. P. (1982). "Reliability, resiliency, and vulnerability criteria for water resource system performance evaluation." Water Resources Research, 18(1), 14-20. Hasofer, A. M., and Lind, N. S. (1974). "Exact invariant second-moment code format." Journal of Engineering Mechanics of the American Society of Civil Engineers, 100( 1), 111-121. Holland, J. H. (1975). Adaptation in natural and artificial systems, MIT Press, Cambridge, MA. Holling, C. S. (1973). "Resilience and Stability of Ecological Systems." Annu. Rev. Ecol. Systems, 4, 1-23. Huang, C. L., and Mayer, A. S. (1997). "Pump-and-treat optimization using well locations and pumping rates as decision variables." Water Resources Research, 33(5), 1001-1012. Jang, Y. S., Sitar, N., and Der Kiureghian, A. (1994). "Reliability analysis of contaminant transport in saturated porous media." Water Resources Research, 30(8), 2435-2448. Krishnakumar, K. (1989). "Micro-genetic algorithms for stationary and non-stationary function optimization." SPIE: Intelligent Control and Adaptive Systems, 1196. 141 Kumar, V. V., Rao, B. V., and Mujumdar, P. P. (1996). "Optimal operation of a multibasin reservoir system." Sadhana Academy Proceedings in Engineering Sciences, 21(4), 487502. Kundzewicz, Z. W. (1989). "Renewal theory criteria of evaluation of water-resource systems. Reliability and resilience." Nordic Hydrology, 20(4-5), 215-230. Leland, D., Anderson, S., and Sterling, D. J. (1997). "The Willamette - a river in peril." Journal of the American Water Works Association, 89(11), 73-83. Lence, B. J., Eheart, J. W., and Brill Jr., E. D. (1990). "Risk equivalent seasonal discharge programs for multidischarger streams." Journal of Water Resources Planning and Management, 116(2), 170-186. Liebman, J. C , and Lynn, W. R. (1966). "The optimal allocation of stream dissolved oxygen." Water Resources Research, 2(3), 581-591. Logan, L. "Performance evaluation for water quality monitored streams." International Symposium on Water Resource Systems Application, Winnipeg, Canada, 133-144. Lohani, B. N., and Saleemi, A. R. (1982). "Recent developments on stochastic programming model for water quality management in rivers." Water Supply and Management, 6(6), 511-520. Lohani, B. N., and Thanh, N. C. (1978). "Stochastic programming model for water quality management in a river." Journal of Water Pollution Control Federation, 50(9), 21752182. Loucks, D. P., Revelle, C. S., and Lynn, W. R. (1967). "Linear programming for water pollution control." Management Science, 14(4), B166-B181. Madsen, H. O., Krenk, S., and Lind, N. C. (1986). Methods of structural safety, Prentice-Hall, Englewood Cliffs, N.J. Maier, H., and Dandy, G. (2000). "Neural networks for the prediction and forecasting of water resources variables: a review of modelling issues and applications." Environmental Modelling and Software, 15, 101-124. Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (1999). "First-order reliability method for estimating reliability, resilience and vulnerability." Water Resources Research, In Press. McKinney, D. C , and Lin, M. D. (1994). "Genetic algorithm solution of groundwatermanagement models." Water Resources Research, 30(6), 1897-1906. 142 Melching, C. S. (1992). "An improved first-order reliability approach for assessing uncertainties in hydrologic modeling." Journal of Hydrology, 132, 157-177. Melching, C. S., and Anmangandla, S. (1992). "Improved first-order uncertainty method for water-quality modeling." Journal of Environmental Engineering, 118(5), 791-805. Melching, C. S., and Flores, H. E. (1999). "Reaeration equations derived from U.S. Geological Survey database." Journal of Environmental Engineering, 125(5), 407-413. Melching, C. S., Yen, C. Y., and Wenzel, H. G. (1990). "A reliability estimation in modeling watershed runoff with uncertainties." Water Resources Research, 26(10), 2275-2286. Melching, C. S., and Yoon, C. G. (1996). "Key sources of uncertainty in QUAL2E model of Passaic River." Journal of Water Resources Planning and Management, 122(2), 105-113. Morgan, D. R., Eheart, J. W., and Valocchi, A. J. (1993). "Aquifer remediation design under uncertainty using a new chance constrained programming technique." Water Resources Research, 29(3), 551-561. Morgan, M. G., and Henrion, M. (1990). Uncertainty: A guide to dealing with uncertainty in quantitative risk and policy analysis, Cambridge University Press, Cambridge. Morshed, J., and Kaluarachchi, J. J. (2000). "Enhancements to genetic algorithm for optimal ground-water management." Journal of Hydrologic Engineering, 5(1), 67-73. Moy, W.-S., Cohon, J. L., and ReVelle, C. S. (1986). "A programming model for analysis of the reliability, resilience, and vulnerability of a water supply reservoir." Water Resources Research, 22(4), 489-498. Mulligan, A. E., and Brown, L. C. (1998). "Genetic algorithms for calibrating water quality models." Journal of Environmental Engineering, 124(3), 202-211. Narayanan, S., and Azarm, S. (1999). "On improving multiobjective genetic algorithms for design optimization." Structural Optimization, 18, 146-155. O'Connor, D. J., and Dobbins, W. E. (1958). "Mechanisms of reaeration in natural streams." Transactions of the American Society of Civil Engineers, 123, 641-684. ODEQ. (1995). "Final issue paper: dissolved oxygen - 1992-1994 water quality standards review." , Oregon Department of Environmental Quality, Portland, Oregon. Oliveira, R., and Loucks, D. P. (1997). "Operating rules for multireservoir systems." Water Resources Research, 33(4), 839-852. Palli, N., Azarm, S., McCluskey, P., and Sundararajan, R. (1998). "An interactive multistage Einequality constraint method for multiple objectives decision making." Journal of Mechanical Design, 120, 678-686. Rackwitz, R. (1976). "Practical probabilistic approaches to design." , Comite European du Beton, Paris, France. Revelle, C. S., Loucks, D. P., and Lynn, W. R. (1968). "Linear programming applied to water quality management." Water Resources Research, 4(1), 1-9. Ritzel, B. J., Eheart, J. W., and Ranjithan, S. (1994). "Using genetic algorithms to solve a multiple-objective groundwater pollution containment-problem." Water Resources Research, 30(5), 1589-1603. Shih, C. S. (1970). "System optimization for river basin water quality management." Journal of the Water Pollution Control Federation, 42(10), 1792-1804. Shinozuka, M. (1983). "Basic analysis of structural safety." Journal of Structural Engineering, 109, 721-740. Sitar, N., Cawlfield, J. D., and Der Kiureghian, A. (1987). "First-order reliability approach to stochastic analysis of subsurface flow and contaminant transport." Water Resources Research, 23(5), 794-804. Skaggs, T. H., and Barry, D. A. (1997). "The first-order reliability method of predicting cumulative mass flux in heterogeneous porous formations." Water Resources Research, 33(6), 1485-1494. Solanki, R. S., Appino, P. A., and Cohon, J. L. (1993). "Approximating the noninferior set in multiobjective linear programming problems." European Journal of Operational Research, 68, 356-373. Solomatine, D. P., Dibike, Y. B., and Kukuric, N. (1999). "Automatic calibration of groundwater models using global optimization techniques." Hydrological Sciences, 44(6), 879-894. Somlyody, L. (1997). "Use of optimization models in water quality planning." Water Science and Technology, 36(5), 209-218. Takyi, A. K. (1991). "Uncertainty analysis for environmental management models: Application of the generalized sensitivity analysis," M.Sc, University of Manitoba, Winnepeg. Takyi, A. K., and Lence, B. J. (1999). "Surface water quality management using a multiplerealization chance constraint method." Water Resources Research, 35(5), 1657-1670. Tetra Tech. (1993). "Willamette River basin water quality study: Willamette River dissolved oxygen modeling component." TC 8983-03, Tetra Tech, Inc., Redmond, WA. Tetra Tech. (1995a). "Willamette River basin water quality study, phase II. Steady-state model refinement component: QUAL2E-UNCAS dissolved oxygen model calibration and verification - Draft Appendix." , Tetra Tech, Inc., Redmond, Washington. Tetra Tech. (1995b). "Willamette River basin water quality study: summary and synthesis of study findings." , Tetra Tech, Inc., Bellevue, WA. Thomann, R. V. (1998). "The future "Golden Age" of predictive models for surface water quality and ecosystem management." Journal of Environmental Engineering, 124(2), 94-103. Tolson, B. A., and Lence, B. J. "Genetic algorithms for reliability-based optimization in water quality management." Annual Conference of The Canadian Society for Civil Engineering, London, Ontario, Canada. Tung, Y.-K. (1990). "Evaluating the probability of violating dissolved oxygen standard." Ecological Modelling, 51, 193-204. Tung, Y.-K., and Hathhorn, W. E. (1988). "Assessment of probability distribution of dissolved oxygen deficit." Journal of Environmental Engineering, 114(6), 1421-1435. Van Note, R. FL, Hebert, R. M., Patel, R. M., Chupek, C , and Feldman, L. (1975). "A guide to the selection of cost-effective wastewater treatment systems." EPA-430/9-75-002, Office of Water Programs Operations, U.S. Environmental Protection Agency, Washington, D.C. Vasquez, J. A., Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (2000). "Achieving water quality system reliability using genetic algorithms." Journal of Environmental Engineering, 126(10), 954-962. Wagner, B. J., and Gorelick, S. M. (1989). "Reliable aquifer remediation in the presence of spatially variable hydraulic conductivity: From data to design." Water Resources Research, 25(10), 2211-2226. Wang, M., and Zheng, C. (1997). "Optimal remediation policy selection under general conditions." Ground Water, 35(5), 757-764. Wang, Q. J. (1997). "Using genetic algorithms to optimise model parameters." Environmental Modelling and Software, 12(1), 27-34. Wardlaw, R., and Sharif, M. (1999). "Evaluation of genetic algorithms for optimal reservoir system operation." Journal of Water Resources Planning and Management, 125(1), 2533. Xu, C , and Goulter, I. C. (1999). "Reliability-based optimal design of water distribution networks." Journal of Water Resources Planning and Management, 125(6), 352-362. Yang, G., Reinstein, L. E., Pai, S., Xu, Z., and Carroll, D. L. (1998). "A new genetic algorithm technique in optimization of permanent 125T p state implants." Medical Physics, 25(12), r0 2308-2315. 145 Zitzler, E., and Thiele, L. (1999). "Multiobjective evolutionary algorithms: a comparative case study and Strength Pareto approach." IEEE Transactions on Evolutionary Computing, 3(4), 257-271. Zongxue, X., Jinno, K., Kawamura, A., Takesaki, S., and Ito, K. (1998). "Performance risk analysis for Fukuoka water supply system." Water Resources Management, 12(1), 13-30. 146 APPENDIX I. Stationarity Plots for the Data used in the Resilience Estimates 180.0 (a) o.o 3 4 5 13-day period (period 1 begins o n J u l y 1) (b) Figure 1-1. Stationarity plot for the (a) mean and (b) standard deviation of the 13-day independent average flows over the low flow season 147 (a) 2.0 Ui o.4 . 0.2 0.0 J , 1 . 2 , 3 , 4 . 5 , 6 7 13-day period (period 1 begins on July 1) (b) Figure 1-2. Stationarity plot for the (a) mean and (b) standard deviation of the 13-day independent average temperature over the low flow season 148 APPENDIX II. Waste Treatment Cost Data Table II-1. Waste treatment cost data for the W W T P s and the P P M s City of Canby City of Wilsonville Cost in Million $US/Year Discrete treatment option 1 2 3 4 5 6 7 8 9 % BOD 5 Removal 80.6 Capital 0.242 5 O&M 0.333 0.575 Total 90.6 0.246 0.347 92.8 94.2 95.4 96.4 0.593 0.338 0.329 0.334 0.405 0.393 0.422 0.731 0.751 0.787 97.1 0.395 98.1 99.1 0.453 0.456 0.453 0.470 0.499 0.697 0.768 Cost in Million $US/Year Discrete treatment % B O D option Removal 1 80.6 2 90.6 3 0.875 0.894 4 5 6 7 1.150 1.223 8 9 City of Corvalis 92.8 94.2 95.4 96.4 O&M Total 0.720 1.421 0.788 1.505 0.913 1.110 1.169 2 90.6 3 4 94.2 0.890 95.4 97.1 0.935 1.101 1.532 1.483 1.462 1.762 1.803 2.046 2.270 2.995 3.244 City of Newberg 0.538 97.1 0.519 0.624 0.638 98.1 99.1 3 4 5 94.2 95.4 97.1 6 7 98.1 99.1 % BOD 1 1.126 1.454 1.590 Capital 0.622 0.635 O&M Total 0.652 0.707 1.273 1.342 0.789 0.821 0.828 0.975 1.322 0.990 1.044 1.610 1.819 2.020 1.316 1.575 2.638 2.865 1.290 Cost in Million $US/Year Discrete Removal 0.830 0.951 0.916 0.923 1.036 1.119 Clackamas Co. Service District #1 Discrete option 0.581 0.607 0.745 Cost in Million $US/Year Discrete treatment % B O D option Removal 1 80.6 2 90.6 Cost in Million $US/Year treatment 0.416 0.472 0.498 0.583 6 Capital 0.702 0.717 98.1 99.1 0.329 0.444 Total 0.717 City of Albany 5 5 6 7 O&M 0.395 0.425 0.452 Cost in Million $US/Year Discrete treatment % B O D option Removal 1 80.6 Capital 0.322 treatment % BOD O&M Total option 80.6 Capital 0.313 Removal 0.388 Capital O&M Total 2 0.701 90.6 1 0.319 80.6 0.409 0.770 0.728 0.778 1.548 3 2 0.432 0.464 0.786 0.857 5 95.4 0.414 0.441 0.490 0.574 0.896 0.904 1.643 4 92.8 94.2 90.6 3 4 94.2 95.4 0.972 1.022 1.015 0.991 1.212 5 1.963 2.234 6 7 96.4 97.1 97.1 0.523 1.204 0.569 1.272 2.476 0.505 6 7 97.7 99.1 1.254 0.595 1.093 1.100 1.493 1.916 2.746 3.558 8 98.1 0.605 9 99.1 0.620 0.816 0.935 5 5 1.643 1.421 1.556 Note: O&M refer to operation and maintenance. 149 Table II-l. Waste treatment cost data for the W W T P s and the P P M s (concluded) City of Salem City of Eugene Cost in Million $US/Year Discrete treatment % B O D option Removal 5 80.6 Capital 2.102 3 90.6 94.2 2.136 2.631 4 5 97.1 97.7 6 99.1 3.215 3.541 5.183 1 2 O&M 2.069 2.395 2.743 3.474 4.620 5.756 Cost in Million $US/Year Discrete treatment % B O D option Removal 1 80.6 5 Total 4.171 4.531 5.374 6.689 8.162 10.938 Tri-City Service District 2 3 4 90.6 94.2 95.4 5 97.1 6 7 98.1 99.1 O&M 0.707 90.6 1.862 3 4 94.2 97.1 2.293 2.050 2.349 5 6 97.7 99.1 2.808 3.098 4.473 3.004 3.959 4.932 3.912 4.642 5.812 7.056 9.405 Discrete treatment % B O D option Removal 1 80.6 Cost in Million $US/Year 5 Total 0.688 0.703 0.872 0.773 0.896 0.918 1.395 1.476 1.767 1.088 1.077 1.491 2.006 1.145 1.432 1.726 2.223 2.924 3.173 1.448 Total 3.612 Oak Lodge Sanitary District 5 Capital O&M 1.778 2 Cost in Million $US/Year Discrete treatment % B O D option Removal 1 80.6 Capital 1.834 2 3 90.6 94.2 4 5 95.4 97.1 6 7 98.1 99.1 Capital 0.412 0.420 O&M 0.471 0.501 0.591 Total 0.883 0.921 0.533 0.553 0.681 0.656 0.834 0.817 1.235 0.735 0.969 1.110 1.390 1.803 1.927 1.125 City of Portland Cost in Million $US/Year Discrete treatment % BOD5 option Removal 1 2 3 4 80.6 90.6 94.2 97.1 Capital 3.587 O&M Total 3.679 3.633 4.484 7.266 4.323 7.956 4.990 5.427 9.474 6.163 11.589 13.820 19.201 5 97.7 6.071 6 99.1 7.749 9.359 9.842 Pulp and Paper Mill Pope and Talbot Inc. James River Paper Co. Inc. Willamette Industries Inc. Smufit Newsprint Corp. (Newberg) Smurfit Newsprint Corp. (Oregon City) Simpson Paper C o . Total Costs in Million $US/Year) Level 1 Level 2 Level 3 Level 4 1.458 1.921 8.589 20.539 0.335 0.374 2.540 5.601 1.340 1.911 5.082 10.915 1.134 1.284 4.376 9.562 _ 0.439 1.894 3.757 0.507 0.565 3.770 8.675 150 APPENDIX III. Correlation Matrices used in the Water Quality Management Models a Table III-l. Original correlation coefficients calculated from the flow and temperature data used for reliability and vulnerability calculation in the water quality management models Variable # Variable # 1 Description 1 2 3 Flow in Willamette River at confluence of Coast and 1 -0.31 0.04 Middle Fork 2 Flow in McKenzie River -0.31 3 Flow in Santiam River 0.04 0.39 1 0.39 1 Table III-2. Original correlation coefficients calculated from the flow and temperature data used for resilience calculation in the water quality management models Variable # Description 1 Flow in Willamette River (t) at confluence of Coast and 1 1 2 Variable # 3 6 7 8 0.14 0.62 0.78 0.15 0.67 Middle Fork 2 Flow in McKenzie River (t) 3 Flow in Santiam River (t) 0.14 1 0.62 0.42 0.42 0.06 0.77 0.26 1 0.24 0.35 0.69 6 Flow in Willamette River (t+1) 0.78 0.06 0.24 7 Flow in McKenzie River (t+1) 0.15 0.77 0.35 0.14 8 Flow in Santiam River (t+1) 1 0.14 0.62 1 0.67 0.26 0.69 0.62 0.42 0.42 1 151
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Genetic algorithms for multi-objective optimization...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Genetic algorithms for multi-objective optimization in water quality management under uncertainty Tolson, Bryan Antony 2000
pdf
Page Metadata
Item Metadata
Title | Genetic algorithms for multi-objective optimization in water quality management under uncertainty |
Creator |
Tolson, Bryan Antony |
Date Issued | 2000 |
Description | This thesis demonstrates the combined usage of a number of novel approaches and techniques in the multi-objective management of water quality systems under uncertainty. The First-Order Reliability Method (FORM) is used to estimate the risk-based system performance indicators of reliability, vulnerability, and resilience that provide measures of the frequency, magnitude and duration of the failure of water resource systems, respectively. FORM accuracy and efficiency for performance indicator estimation is compared extensively with Monte Carlo Simulation (MCS). Genetic Algorithms (GAs) are demonstrated as a robust optimization technique by solving various multi-objective water quality management models that optimize the performance indicators and the total point source waste treatment cost. In addition, the Tradeoff Surface Representation (TSR) Algorithm is incorporated as a general multi-objective technique for accurate and efficient identification of convex tradeoff surfaces. The Willamette River Basin in Oregon, USA is utilized as the water quality management case study for the demonstration of all techniques. The performance indicators are estimated with respect to meeting dissolved oxygen (DO) standards and ambient DO is simulated using a QUAL2E water quality response model. Results show that FORM estimates of the performance indicators, while significantly less accurate than MCS estimates, seem to provide reasonable results when utilized within the multiobjective water quality management models. A comparison of FORM and MCS shows that while FORM is more efficient relative to MCS, the difference in efficiencies is significantly less than previously reported in the literature. The TSR Algorithm, in comparison with the commonly used Constraint Method for multi-objective tradeoff curve generation, is shown to produce a superior representation of the tradeoff curve. Furthermore, the TSR Algorithm is also shown to produce a maximum amount of auxiliary information regarding the bounds on the location of the tradeoff curve between tradeoff points. |
Extent | 8850703 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-07-20 |
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.0064020 |
URI | http://hdl.handle.net/2429/10987 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2000-0596.pdf [ 8.44MB ]
- Metadata
- JSON: 831-1.0064020.json
- JSON-LD: 831-1.0064020-ld.json
- RDF/XML (Pretty): 831-1.0064020-rdf.xml
- RDF/JSON: 831-1.0064020-rdf.json
- Turtle: 831-1.0064020-turtle.txt
- N-Triples: 831-1.0064020-rdf-ntriples.txt
- Original Record: 831-1.0064020-source.json
- Full Text
- 831-1.0064020-fulltext.txt
- Citation
- 831-1.0064020.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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064020/manifest