ESTIMATING RESILIENCE FOR HYDROTECHNICAL SYSTEMS by YiLi BEng, Chongqing Jiaotong University, 1994 MSc, University of Guelph, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENT FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA October 2007 © Yi Li, 2007 ABSTRACT With the cycle of natural conditions (hydrological or meteorological), a repairable hydrotechnical system may have a recovery capacity after a system failure. Resilience is a probabilistic performance indicator that measures such a recovery capacity. It is also one of performance indicators that can be used to measure the sustainability of a water resources project. In nature, resilience is a conditional probability that, given a system failure at an earlier time moment, the system can resume functioning at a later time moment. The conditional probability is a close derivative of two-state transition probabilities in a stochastic process. The current estimating techniques for the two-state lag-1 transition probabilities may require a high computational cost, or a complex implementation procedure. This thesis explores practical approximation methods for estimating the two-state lag-1 transition probabilities in discrete processes, and the resultant resilience. Two methods are proposed herein. By developing a concise and straightforward approximation for the lag-1 autocorrelation coefficient of system performance functions at two consecutive moments on the basis of the first-order reliability method, one method improves the conventional approaches that estimate the transition probabilities based on a bivariate normal distribution of the system performance functions. The other employs a linear stochastic prediction of the performance function based on the mean point of failure or safety domain. The mean points are one of the fundamental properties of reliability problem. The thesis mathematically defines the mean points, and develops a practical approach for approximating them based on the first-order reliability method. The mean points-based method improves the conventional approaches that estimate the transition probabilities based on a linear stochastic prediction of performance function. The two methods developed herein are demonstrated in a typical river hydrology example for estimating seasonal lag-1 resilience. The proposed approximation methods for estimating the transition probabilities in discrete processes can be also applicable to other civil engineering problems, such as marine/offshore steel structure fatigue analysis, river navigability study, or concrete reinforcement corrosion prediction. n TABLE OF CONTENTS Abstract ii Table of contents iii Table of tables : vi Table of figures vii Preface viii Acknowledgements ix Co-authorship statement x 1 Introduction 1 1.1 Background 2 1.2 Resilience ., 5 1.3 Thesis objectives 6 1.4 Research methods 8 1.5 Structure of thesis 10 References 12 2 Mean points of failure and safety domains: estimation and application 15 Preface 16 2.1 Introduction 17 2.2 Definitions 17 2.3 First-order mean points of failure and safety domains 18 2.3.1 FORM 18 2.3.2 First-order approximations of MPSD and MPFD 19 2.3.2.1 MPFD 20 2.3.2.2 MPSD 21 2.3.3 Numerical examples 22 2.4 Estimates of two-state transition probabilities ; 24 2.4.1 Two-state transition probabilities 24 2.4.2 A method of estimating transition probabilities based on MPSD and MPFD 27 2.4.2.1 A stationary linear univariate example 27 2.4.2.2 A stationary nonlinear univariate example 32 2.4.3 Summary : 34 2.5 Conclusions 35 References 36 3 Estimating resilience for water resources systems 38 Preface 39 3.1 Introduction.. 40 3.2 Background 41 3.2.1 Definitions 41 3.2.2 Previous work 42 3.3 Proposed methodologies 44 3.3.1 The bivariate normal distribution (BND) method 45 3.3.2 The stochastic linear prediction (SLP) method 48 3.3.3 Summary 51 3.4 Numerical example 52 3.4.1 Problem setting: water supply from two sources 52 3.4.2 Estimates of reliability and resilience 54 3.4.3 Results 55 3.5 Conclusions 58 References : 60 4 On the applicability of Rice's formula in stochastic hydrological modeling 63 Preface 64 4.1 Introduction 65 4.2 Rice's formula 65 4.3 Crossings in discrete processes 66 4.4 Numerical examples 67 4.5 Conclusions 70 References 71 5 Conclusions . 72 5.1 Summary 73 5.2 Conclusions 76 5.3 Future work 77 iv 5.4 Limitations 78 References 79 A p p e n d i x A N u m e r i c a l a p p r o x i m a t i o n s o f d e s i g n p o i n t s i n r e l i a b i l i t y a n a l y s i s u n d e r p a r a m e t r i c c h a n g e s 8 0 Preface 81 A. 1 Introduction 82 A.2 Parametric sensitivity measures 84 A. 3 Numerical approximation algorithms 85 A.3.1 Numerical schemes 86 A.3.2 Computational procedures 88 A.4 Numerical examples 91 A.4.1 Example 1: A plane frame structure 92 A.4.2 Example 2: A reinforced concrete beam 94 A.4.3 Example 3: A rectangular column 95 A.4.4 Example 4: Riverine dissolved oxygen response 96 A.4.5 Example 5: Impacts of noise 97 A.5 Conclusions 98 References 100 V LIST OF TABLES Table 2.1. Statistical properties of random variables in numerical examples 22 Table 2.2. Computational results of PFS in the nonlinear example 34 Table 3.1. Monthly water demands and mean river flows (Mm3) 52 Table 3.2. Estimates of Ps and Re\ 56 Table 3.3. Number of performance function evaluations 56 Table 3.4. Impacts of an additional water supply 58 Table A.l . Statistical properties of random variables in numerical examples 93 Table A.2. Computational results for Example 1: A plane frame structure 93 Table A.3. Computational results for Example 2: A reinforced concrete beam 95 Table A.4. Computational results for Example 3: A rectangular column 95 Table A.5. Computational results for Example 4: Riverine dissolved oxygen response 96 Table A.6. Computational results for Example 5: Impacts of noise 97 vi L I S T O F F I G U R E S Figure 2.1. Schematic description of FORM in a bivariate u-space 19 Figure 2.2. Illustrations of z and its normal distribution 20 Figure 2.3. ZF and zs as function of zp 21 Figure 2.4. MPFD and MPSD in numerical example 23 Figure 2.5. MPFD and MPSD in numerical example 24 Figure 2.6. A stationary Gaussian process V 28 Figure 2.7. Sample paths of Vbetween to and to+l 29 Figure 2.8. Computational experiments on PFS(*O) (C = 0) 30 Figure 2.9. Computational experiments on PFSOO) (C = -0.5, -1.0, -1.5) 30 Figure 2.10. Computational experiments on Re(t0) (C = 0 to -3.0) 31 Figure 3.1. A probabilistic space of z(ti) and z(ti) 46 Figure 3.2. A bivariate example of FORM at two consecutive moments in time 49 Figure 3.3. Sample paths of u inherent in Pfs(fi, ti) and Re\(t\) 51 Figure 3.4. MCS estimates and their 95% confidence intervals 56 Figure 4.1. Examples of Rxif) 66 Figure 4.2. A stationary Gaussian process X 68 Figure 4.3. Computational results for Pu\(to) 69 Figure 4.4. Computational results fori?ei(fo) 69 Figure A.l . Configuration of initial-value problem in a bivariate u-space 86 Figure A.2. The working flowchart of J-Euler and J-iEuler 89 Figure A.3. The operation procedure of S-iEuler 90 Figure A.4. Illustration of efficiency comparison between HL-RF and numerical techniques.... 90 Figure A.5. A plane frame structure 92 PREFACE This thesis has led to the following publications: • Li, Y., and Lence, B. J. (2007b). "Mean Points of Failure and Safety Domains: Estimation and Application." Probabilistic Engineering Mechanics, 22(3), 257-266. • Li, Y., and Lence, B. J. (2007a). "Estimating Resilience for Water Resources Systems." Water Resources Research, 43(7), W07422. Li, Y., and Lence, B. J. (2007c). "Numerical Approximations of Design Points in Reliability Analysis under Parametric Changes." ASCE Journal of Engineering Mechanics, in press. • Li, Y., and Lence, B. J. (2006). "On the Applicability of Rice's Formula in Stochastic Hydrological Modeling." ASCE Journal of Hydrological Engineering, in review. • Li, Y., and Lence, B. J. (2005). "On Risk Analysis of Water Resources Systems under Non-Stationary Conditions." Proceedings of World Water & Environmental Resources Congress 2005, ASCE, Anchorage, Alaska. In addition, another conference paper was produced during the author's PhD study, though it is not directly related to this thesis: • Li, Y., Tang, Y., and Li, D. (2006). "Hydroelectric potential of stormwater runoff from the University of British Columbia's Vancouver campus." BCWWA Annual Conference and CAWQ Western Regional Symposium, Whistler, British Columbia, Canada. v i a A C K N O W L E D G M E N T S Here I want to express my appreciation to all of the people who have been supporting and helping me through my PhD study. I would like to thank my supervisor, Dr. Barbara Lence. Her constant supports, kind guidance, and informative instructions are a key basis for this thesis, always bringing me inspirations. Especially, my writing skills are substantially improved under her trainings. Strong supports from my family are my spiritual prop. I highly appreciate my wife Qiang for her great understanding during my graduate studies. She has been making huge efforts to look after our lovely daughter Cheng Cheng. She even partly sacrificed her career in order to support me. Also, Cheng Cheng sacrifices much happy playing time with father. Qiang and Cheng Cheng, I love you two. I owe too much to them, as well as my parents and parents-in-law. Thanks are also given to my comprehensive examination committee, advisory committee, and final examination committee, including Dr. Younes Alila, Dr. Sue Baldwin, Dr. Scott Dunbar, Dr. Ricardo Foschi, Dr. Terje Haukaas, Dr. Bernard Laval, Dr. Edward McBean, Dr. Robert Millar, and Dr. Siefried Stiemer. This thesis benefits from their constructive comments and suggestions. Also, all of my friends here provide appreciable help to me. I cherish those unforgettable hours I spent in working with them and sharing joy with them. Especially, I am grateful to our research mates, Ali, Andrew, and Hamid, for informative discussion. My personal effort is another key to the success of this research. The careful and creative work gives me confidence. I appreciate all financial supports for my PhD study, NSERC PGS D, NSERC Discovery Research Grant, PARC Supplement to NSERC Scholarships, UBC Graduate Entrance Scholarship and PhD Tuition Fee Award, UBC Bursary, Earl R. Peterson Memorial Scholarship in Civil Engineering, and Canada Study Grant for Students with Dependents. These supports not only financially aid me, but also, more importantly, spiritually encourage me to complete my study. IX CO-AUTHORSHIP S T A T E M E N T Yi Li was the lead and principal researcher of the work contained in the thesis titled "Estimating Resilience for Hydrotechnical Systems". Dr. Barbara Lence, the Research Supervisor of the thesis, provided inspiration and supported the writing of the papers included in the thesis. CHAPTER 1 INTRODUCTION 1.1 BACKGROUND Uncertainty is one of the primary challenges faced by engineers in the planning, design, construction, operation, and maintenance of major civil engineering projects. The following are several examples. In planning a municipal stormwater system, the rainfall intensity is a key design parameter. Since it is highly uncertain, intensity-duration-frequency curves are widely used in order to identify an appropriate intensity for a given urban watershed. For breakwater designs, wave height plays an important role in sizing armor units and in assessing wave runup and wave overtopping. To account for the natural variation of local wave height, wave statistics are developed in practice, and used to identify an appropriate design wave height for a specific breakwater structure. In structural designs, the engineer needs the important information on material strength, as well as loads and their combinations, but these are never exactly known in reality. To mitigate the effects of these uncertainties, most contemporary structural codes impose reliability-based factors on material strength and loads. These are a rather limited number of examples of uncertainty, and there are many more in engineering practice. Civil engineering design is widely founded in the process of making decision under various uncertainties. Clearly, uncertainty analysis may significantly aid in this process. According to Ayyub and Chao (1998), uncertainties in engineering systems may fundamentally stem from ambiguity and vagueness. Ambiguity-related uncertainties mainly result from non-cognitive sources, such as physical randomness, statistical errors, lack of knowledge, and modeling errors. Vagueness-related uncertainties mainly result from cognitive sources, such as the inappropriate definition of a problem, human factors, and the inappropriate definition of relationships between parameters in a problem. In practice, ambiguity-related uncertainties are generally analyzed with probability theory, while vagueness-related uncertainties are generally analyzed with fuzzy set theory. This thesis focuses on the ambiguity-related uncertainties in hydrotechnical systems. More specifically, probabilistic uncertainty analyses for hydrotechnical systems with random system inputs are undertaken. In a system with a series of random inputs, the system output is also a random variable. The central task of uncertainty analysis is to estimate the statistical properties of the system output. Usually, the analytical solution of the task is not available, and practitioners are interested in approximate solutions. Depending on the problem, practitioners generally analyze uncertainties for a system in two different but related ways, i.e. by 1) developing the probabilistic distribution 2 or statistical moments of system output, and 2) evaluating performance indicators (Pis) of system output. A very typical scenario of uncertainty analysis is the derivation of the probabilistic distribution of the system output, i.e., the probability density function (PDF) and the cumulative density function (CDF). For instance, the river flow may be a predominate uncertain factor in the design of a local dyke system. To implement an uncertainty analysis for the river flow, a regular hydrological practice is to fit a probabilistic distribution to the available annual maximum flowrate data series, in other words a PDF/CDF pattern. However, both functions are usually difficult, if not impossible, to derive due to the complexity of real problems, such as multiple random variables and nonlinearity. In these cases, statistical moments of the system output, particularly second-order moments (mean and standard deviation) may be the major interest of the decision maker. Many uncertainty analyses in practice are devoted to simply the derivation of the second-order moments of system output. As an example, such an analysis helps to identify important uncertain sources in the riverine dissolved oxygen (DO) prediction (Chadderton et al. 1982; Tung and Hathhorn 1988; Melching and Yoon 1996). For further information on developing the second-order moments of system output, see Tung (1996), Tyagi and Haan (2001), and Haan (2002). Oftentimes, the decision maker may want to see the likely behavior of a system in terms of a specific aspect under random and varied deterministic input conditions. The likely system behavior may be described by various Pis, which are developed within a stochastic process framework. Many Pis have been developed to reflect different aspects of system performance, depending on the "feature" of a system relative to its failure mechanism. In routine civil engineering systems, such a feature is perhaps often related to their natural repairablity after failure. Characterized by structural failure, a naturally non-repairable system ceases to work once failed. In terms of a dam breach or a levee slide, as an example, the function of the dam or the levee is considered to terminate permanently. Characterized by performance failure, naturally repairable systems are able to naturally recover from failures through the ability of the system to assimilate stress or the cycle of stochastic inputs; e.g., most water resources systems featured by meteorological and hydrological cycles recover naturally from extreme events. For example, a regional drought occurs due to poor meteorological conditions, but will be eliminated by following precipitation. Water quality deterioration caused by extreme summer low flows will be alleviated in the following rainy season. A river flood, which may be catastrophic, will 3 eventually recede. For further information regarding the use of Pis in repairable and non-repairable hydrotechnical systems, see Li and Lence (2005). The most common PI is the availability (or instantaneous reliability) function that measures the probability of the success (or safe, safety, satisfactory) system state at any single moment. (The opposite system state is called the failure or unsatisfactory state.) This PI is more useful in repairable systems. In reservoir regulation, the likelihood of adequate water supply in critical months, i.e., the instantaneous reliability, is a major consideration in sizing reservoir and in designing operational policies. Different policies may apply for dry or wet months. The availability function is the foundation of all other Pis. It should be noted that the value of the availability function at a single moment is customarily also called reliability. This terminological convention should not be confused with a related but different PI, the reliability function. For further details on reliability methods, see relevant textbooks (Ang and Tang 1984; Madsen et al. 1986; Sundararajan 1995; Tung et al. 2006). The reliability function measures the probability of success of a system during a specific period. It in fact demonstrates the "cumulative" risk for a system. The reliability function is a PI useful in both non-repairable and repairable systems. In repairable hydrotechnical systems, the reliability function is useful in evaluating the flooding risk for a specific short- or long-term period (Tung and Mays 1980; Gui et al. 1998; Yanmaz 2000). The reliability function is especially important to non-repairable systems, for which the failure will be permanent. It is useful in evaluating time-dependent structural reliability in terms of aging material (Andrieu-Renaud et al. 2004; Li et al. 2005), dynamic loads such as earthquake and wind (Der Kiureghian 2000), or service period requirement (USACE 2006). A number of other Pis have also been developed for measuring transition probabilities in the time domain between different system states. The transition probabilities are useful in both non-repairable and repairable systems. A reasonable maintenance program of infrastructures may be scheduled on the basis of transition probabilities between various physical condition states of infrastructures during their deterioration process (Madanat and Ibrahim 1995; DeStefano and Grivas 1998; Mishalani and Madanat 2002), and transition probabilities between the flooding and the normal conditions, or between the drought and the normal conditions, at a given time moment are one of key topics in stochastic hydrology. The transition probabilities also play an important role in evaluating the reliability function (Kapur and Lamberson 1977), This thesis focuses on the one-step (or lag-1) transition probabilities between two opposite states in the discrete time domain, as well as one of the derivatives of the transition probabilities, resilience, 4 which can be used to measure the likely duration of a failure in repairable systems. More details regarding resilience are provided in the next section. In summary, the focus of this thesis is on the probabilistic study of lag-1 two-state transition probabilities and related lag-1 resilience of hydrotechnical systems. 1.2 RESILIENCE Originally, the concept of resilience was developed for ecological systems, indicating the ability of a system to absorb various changes due to forcing factors, and then to persist (Holling 1973; Fiering 1976). A related concept is stability, which represents the ability of a system to stabilize around its equilibrium state in the face of various stresses. A high stability usually corresponds to a low resilience. The original physical meaning of resilience implies that, after a system survives, or absorbs, certain changes, its equilibrium state may correspondingly change. The concept was readily adapted for water resources engineering as a measure of the system's recovery ability (Hashimoto et al. 1982). In the event of a failure, a system which is able to recover fast typically experiences less damage, and is thus less vulnerable to external disturbance. In a regional water supply system, for example, a two-day water shortage obviously causes a more serious economical loss than a one-day shortage. An aquatic environment may survive serious water quality deterioration for a short period, but a long exposure to pollutants can cause irreversible damage to the environment. Today, resilience is accepted as one of three important Pis that can be jointly used to measure the sustainability of a water resources project (ASCE Task Committee on Sustainability Criteria 1998; UNESCO/DTP 1999). The other two Pis are reliability and vulnerability, which measure the likelihood of a system success and the likely magnitude of a failure event, respectively. Generally, resilience is defined in practice as the conditional probability that, given a system failure at a given moment, the system will recover at a later moment. In other words, resilience indicates the probability of a transition from failure state to success state between two points in time, consecutive or inconsecutive. The time lag between both moments could be either considerable or infinitesimal. In hydrotechnical practice, the time lag is usually a single time step, leading to the lag-1 resilience. More details on resilience are provided in Chapters 2 and 3. The lag-1 resilience has important implications in hydrotechnical system planning. In a regional water supply system, the lag-1 resilience, along with reliability and vulnerability, helps in establishing effective and reasonable emergency policies and alternatives for water use and supply. For example, assume that the available water in May is high due to spring runoff, but 5 high demand and short supply in the following months result in a low likelihood of being able to meet demand in those months. In this case the reliability in May is high, but the lag-1 resilience for the following months may be considered to be low. Here, the utility is advised to implement stringent water use policies in May as well as the following months, and source additional supplies as needed in June and July to guard against shortages in those months. Sampling methods (or simulations) may be used to estimate the lag-1 resilience of a system through generating synthetic data series with a stochastic model for all random variables involved, and then counting the resultant output series. Sometimes, real historical data records that are considered adequately representative are directly used for describing random variables. Because the lag-1 resilience is a conditional probability, however, simulations may exact a high computational cost, as shown in Chapter 3 and Li and Lence (2007a). Appropriate approximations may considerably reduce computational efforts in estimating resilience. In practice, the main approximation method for estimating lag-1 two-state transition probabilities and related resilience is developed based on a bivariate normal distribution (BND) (Der Kiureghian 2000; Maier et al. 2001; Andrieu-Renaud et al. 2004). However, the current techniques employed in the BND-based method for estimating the correlation coefficient between consecutive performance functions may significantly increase the total number of random variables required for probabilistic evaluation, thus complicating the solution procedure, as discussed in Chapter 2 and Li and Lence (2007b). Another applied method is developed based on a stochastic linear prediction (SLP) of performance function. At present, the SLP is implemented using Rice's formula (Rice 1944, 1945; Melchers 1999; Li et al. 2005). However, Rice's formula-based approaches are not very suitable for discrete stochastic processes that typify water resources systems, as shown in Chapter 4 and Li and Lence (2006). Hence, more practical approximation methods are needed for estimating lag-1 resilience. This thesis is devoted to the exploration of such practical approximation methods. 1.3 T H E S I S O B J E C T I V E S The objectives of this thesis are: To develop practical approximation methods for estimating lag-1 two-state transition probabilities in discrete processes; • To apply the approximation methods for estimating lag-1 resilience; To test the developed resilience methods in a typical hydrotechnical case. 6 The methods developed herein are intended to be capable of accommodating common stochastic conditions in water resources systems, i.e., multiple random variables with non-Gaussian (or normal) distributions, nonlinear performance functions, and non-stationarity. In some water resources problems, the system output is a linear function of a single random variable that follows a normal process. For example, the stream flowrate may be the only variable used in evaluating regional flooding or drought risk (Kundzewicz 1989; Salas et al. 2001; Abi-Zeid et al. 2004; Loaiciga 2005). Thus Pis for such issues are investigated based solely on the flowrate dataset of a desirable stream. Many water resources systems might typically involve multiple influential random variables that follow non-Gaussian processes, and their system performance function could be a nonlinear function of these random variables. For example, the interaction between carbonaceous biochemical oxygen demand and DO in stream is a complicated mechanism, which is a nonlinear function of several principal uncertain sources such as stream flowrate and water temperature. Hence it is necessary to accommodate the multiple variables and non-linearity in the developed resilience methods. Built on individual stochastic processes, the system output is also a stochastic process. A stochastic process is a series of dynamic generalizations of random variables, and the stationarity of a stochastic process becomes important in analyzing Pis. If the system output process is stationary, e.g., its statistical properties are constant over time, the estimated Pis will also remain constant over time, thus the evaluation of Pis is significantly simplified. As an example, the conventional frequency analysis assumes an annual maximum stream flow series to be stationary. On the other hand, in a non-stationary system output process, the statistical properties of random variables vary over time. Consequently, the estimated Pis of the system will also vary over time, thus complicating the evaluation procedure. Non-stationarity is quite typical in hydrotechnical systems. Hydrometric data may experience long-term non-stationarity, such as shifts or annual trends due to artificial watershed regulations (through reservoir, dam, or canalization), urbanization, industrialization, or climate change. Short-tern non-stationarity is typified by various seasonal periodicities due to the annual meteorological cycle, such as seasonal, monthly, or weekly variations of hydrological conditions. It is thus important to consider non-stationary hydrological conditions, especially short-term, in system risk analyses (Deo et al. 1997; Weiler et al. 2000; Durrans et al. 2003; Coles and Tawn 2005). Hence it is necessary for the developed resilience methods to be able to consider the non-stationarity. The resilience methods developed herein may be practically applied in the system planning of water resources projects. A primary application may be the short-term (seasonal, monthly, 7 biweekly, or weekly) operations of a repairable system featured mainly by periodic hydrological conditions. Typical examples include water quality management in lakes and rivers, and integrated water use management. In these systems, the short-term variations of resilience may be vital information in understanding system behaviors, and in evaluating tradeoffs in risk mitigation. A typical water supply example is presented in Chapter 3. This example demonstrates system behaviors using resilience, and investigates the effects of different schemes for the addition of water source to the existing system. Extended applications of the approximation methods for estimating lag-1 two-state transition probabilities are expected. For instance, the methods may be used in structural reliability analysis of systems facing dynamic conditions and in the strength prediction of aging material, where transition probabilities are a key factor in evaluating the reliability function. 1.4 R E S E A R C H M E T H O D S Fundamental principles that applied in this research are outlined here. Uncertainties to be studied As mentioned previously, this research concentrates on the ambiguity-related uncertainties; that is, uncertain variables are assumed to be random, and can be addressed statistically. In water resources engineering analyses, major random uncertainty sources are typically linked to physical randomness, or natural variations such as stream flowrates. The effect of the remaining random uncertainty sources such as modeling and statistical errors may be reduced through appropriate human efforts in practice (Tung 1996). Focusing on the development of practical resilience estimation techniques, this research will study natural variations as random variables. It is worth stressing the necessity and reasonableness of statistical treatments for natural variations in water resources systems risk analyses. Water resources systems are extremely complex due to meteorological and geophysical effects on hydrological phenomena, such as the rainfall-runoff process, transport mechanisms, and evaporation. In practice, many variables are considered to suffer uncertainties merely because humans are unable to adequately understand their inner physical mechanisms in a deterministic manner, not because they are truly random. They are treated as random variables for convenience. Furthermore, stochastic process models, which are pure mathematical mechanisms merely constructed on data, without any physical representation, are introduced to mimic the historical development of those so-called random 8 variables. One of the most common applications of stochastic process models is the well-known flood frequency analysis. Although widely accepted practice in engineering analyses, unfortunately, these statistical treatments may not be theoretically justified (KIernes 2000a, 2000b) Without any other justifiable tools, however, it is necessary to compromise by considering those variables to be purely random, and then by applying stochastic process models to analyze them and their effects on system performance. Therefore, although theoretically unreasonable more or less, these statistical treatments may be practically necessary. Practitioners must recognize that they indeed are a choice due to no choice, though their implementation might be enforced in nearly every project. This research takes these statistical treatments as a fundamental theoretical basis, specifically assuming probabilistic distributions and stochastic process models for random variables. Stochastic models to be used A stochastic process model will be applied to mimic the physical driving mechanism, or the historical development, of random variables. Common models in hydrological practice include point process models and time series models. This research will use the time series technique for stochastic process modeling of random variables. The point process mimics occurrences of an event as discrete, random points over the time horizon, ignoring all details between these points. The event is identified by a specified threshold or standard. Widely applied point process models are those based on the renewal process, which is typified by the well-known Poisson process. The point process is very effective in modeling linear systems with a univariate random variable, i.e., low- or high-flow events identified by a given flowrate. In complex non-linear multivariate systems, however, it may not apply. Also, the construction of a point process is too sensitive; that is, slightly different thresholds may result in different model structures. For hydrological applications of the point process, see North (1980), Lee et al. (1986), Chang (1987; 1989) and Abi-Zeid et al. (2004). The time series modeling technique is capable of reproducing the most important statistical features of many hydrological data time series (Salas 1996). It divides the original data sequence of a random variable into deterministic components and a random component. The former are deterministic functions of time. The latter, also called residual or noise, forms a stationary series, which can be fitted with a probabilistic distribution. Thus the original series can be forecasted through forecasting the residual. The essence of the time series technique is in analyzing the residual. A basic time series model for the residual series is the autoregressive moving-average 9 (ARMA) model for a single variable, a linear time series model. A widely used extension of A R M A in the multivariate situation is vector or multivariate A R M A ( V A R M A ) model, which is built based on the interdependence structure of a set of related processes. Main uses of the technique lie in short-term forecasting for real-time operations, and the generation of long-term synthetic data for planning purposes. Obviously, this technique is more powerful than the point process, especially considering that it can be used to study non-linear combinations of multivariate variables. Moreover, it has been successfully applied in hydrological modeling over four decades, and is a well-guided tool today (Yevjevich 1972; Kottegota 1980; Salas et al. 1980; Brockwell and Davis 1996; Salas 1996). This thesis develops approximation methods for estimating lag-1 two-state transition probabilities in discrete processes based on V A R M A . 1.5 S T R U C T U R E O F T H E S I S This thesis is organized as a series of three manuscripts that have been accepted or are being reviewed for publication in peer-reviewed journals. As a paper manuscript, each of the three chapters (Chapters 2 to 4) can be read and understood independently. Each of them contains a preface that links that chapter to the others, a main body, and references. In each of the chapters, figures and tables are embedded in the main body. Chapter 5 presents conclusions and future work. In addition, Appendix A presents a piece of work that is related to general reliability analysis, but not directly related to resilience. A l l Chapters and Appendix A adopt the reference style applied in American Society of Civi l Engineers journals. Chapter 2 introduces the mean points of the failure and safety domains, based on which an innovative method for estimating lag-1 two-state transition probabilities in discrete processes is developed. The method employs an SLP of the performance function at a later moment based on the mean points of the failure and safety domains at an earlier moment. The method is then tested in systems constructed on univariate normal processes for estimating the transition probabilities and lag-1 resilience. A manuscript resulting from this Chapter has been published in Probabilistic Engineering Mechanics, Elsevier, and referenced as L i and Lence (2007b). In Chapter 3, two new approaches for estimating lag-1 resilience are developed. The first one is a BND-based method, in which a concise but effective approximation for the correlation coefficient of the performance function is formulated. Developed based on the mean points of the failure and safety domains, the second one is a generalized form of the application of the mean points in Chapter 2. The two approaches are demonstrated for a water supply problem. A 10 manuscript resulting from this Chapter has been published in Water Resources Research, American Geophysical Union, and referenced as Li and Lence (2007a). Chapter 4 examines in depth the application of the traditional SLP-based method using Rice's formula for estimating the lag-1 two-state transition probabilities in discrete stochastic process modeling. The results show that the Rice's formula-based method may be an inappropriate modeling tool in hydrological applications. This examination provides insights into the lag-1 two-state transition probabilities and lag-1 resilience in discrete processes. A manuscript resulting from this Chapter is being peer-reviewed in the Journal of Hydrological Engineering, American Society of Civil Engineers, and referenced as Li and Lence (2006). The final chapter, Chapter 5 summarizes conclusions and limitations of the methods developed herein, and discusses future work. Appendix A explores a new strategy for repeatedly estimating reliability under frequent parameter variations. The central idea is to update the design point in the parameter domain, rather than in the traditional random variable domain, by evaluating several parametric sensitivity measures, which are systems of nonlinear first-order ordinary differential equations relating the design point to parameter changes. Appendix A is a byproduct of this thesis, funded by the author's scholarships. A manuscript resulting from this work has been accepted for publication in the Journal of Engineering Mechanics, American Society of Civil Engineers, and referenced as Li and Lence (2007c). 11 REFERENCES Abi-Zeid, I., Parent, E., and Bobee, B. (2004). "The stochastic modeling of low flows by the alternating point processes approach: Methodology and application." Journal of Hydrology, 285,41-61. Andrieu-Renaud, C , Sudret, B., and Lemaire, M. (2004). "The PHI2 method: a way to compute time-variant reliability." ReliabUity Engineering & System Safety, 84(1), 75-86. Ang, A. H. S., and Tang, W. H. (1984). Probability Concepts in Engineering Planning and Design, Vol. II: Decision, Risk, and Reliability, John Wiley, New York. ASCE Task Committee on Sustainability Criteria. (1998). Sustainability Criteria for Water Resource Systems, American Society of Civil Engineers. Ayyub, B. M., and Chao, R.-J. (1998). "Uncertainty Types and Modeling Methods." Uncertainty Modeling and Analysis in Civil Engineering, B. M. Ayyub, ed., CRC Press LLC, Boca Raton, FL, 1-32. Brockwell, P. J., and Davis, R. A. (1996). Introduction to Time Series and Forecasting, Springer-Verlag, New York, NY. Chadderton, R. A., Miller, A. C , and McDonnell, A. J. (1982). "Uncertainty Analysis of Dissolved Oxygen Model." ASCE Journal of the Environmental Engineering Division, 108(EE5), 1003-1013. Chang, T. J. (1987). "Drought Analysis in the Ohio River Basin." ASCE Proceedings of Engineering Hydrology, 601-609. Chang, T. J. (1989). "Characteristics of Streamflow Drought." New Directions for Surface water Modeling, IAHS Baltimore Symposium, IAHS, Baltimore, 333-341. Coles, S., and Tawn, J. (2005). "Seasonal Effects of Extreme Surges." Stochastic Environmental Research and Risk Assessment, 19(6), 417-427. Deo, M. C , Sherief, M., and Sarkar, A. (1997). "Wave Height Estimation Using Disaggregation Models." ASCE Journal of Waterway, Port, Coastal, and Ocean Engineering, 123(2), 63-67. Der Kiureghian, A. (2000). "The Geometry of Random Vibrations and Solutions by FORM and SORM." Probabilistic Engineering Mechanics, 15(1), 81-90. DeStefano, P. D., and Grivas, D. A. (1998). "Method for Estimating Transition Probability in Bridge Deterioration Models." ASCE Journal of Infrastructure Systems, 4(2), 56-62. Durrans, S. R., Eiffe, M. A., Thomas Jr., W. O., and Goranflo, H. M. (2003). "Joint Seasonal/Annual Flood Frequency Analysis." ASCE Journal of Hydrological Engineering, 8(4), 181-189. Fiering, M. B. (1976). "The Role of Systems Analysis in Water Program Development." Natural Resources Journal, 16,759-771. Gui, S., Zhang, R., and Xue, X. (1998). "Overtopping Reliability Models for River Levee." ASCE Journal of Hydraulic Engineering, 124(12), 1227-1234. Haan, C. T. (2002). Statistical Methods in Hydrology, Iowa State Press, Ames. 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. Holling, C. S. (1973). "Resilience and Stability of Ecological Systems." Annual Review of Ecological Systems, 4, 1-23. Kapur, K. C , and Lamberson, L. R. (1977). Reliability in Engineering Design, John Wiley & Sons, Inc. Klemes, V. (2000a). "Tall Tales about Tails of Hydrological Distribution, I." ASCE Journal of Hydrological Engineering, 5(3), 227-231. 12 Klemes, V. (2000b). "Tall Tales about Tails of Hydrological Distribution, II." ASCE Journal of Hydrological Engineering, 5(3), 232-239. Kottegota, N. T. (1980). Stochastic Water Resources Technology, The Macmillan Press Ltd., Hong Kong. Kundzewicz, Z. (1989). "Renewal Theory Criteria of Evaluation of Water-Resources Systems: Reliability and Resilience." Nordic Hydrology, 20, 215-230. Lee, K. S., Sadeghipour, J., and Dracup, J. A. (1986). "An Approach for Frequency Analysis of Multiyear Drought Durations." Water Resources Research, 22(5), 655-662. Li, C. Q., Lawanwisut, W., and Zheng, J. J. (2005). "Time-Dependent Reliability Method to Assess the Serviceability of Corrosion-Affected Concrete Structures." ASCE Journal of Structural Engineering, 131(11), 1674-1680. Li, Y., and Lence, B. J. (2005). "On Risk Analysis of Water Resources Systems under Non-Stationary Conditions." World Water & Environmental Resources Congress 2005, ASCE, Anchorage, Alaska. Li, Y., and Lence, B. J. (2006). "On the Applicability of Rice's Formula in Stochastic Hydrological Modeling." ASCE Journal of Hydrological Engineering, in review. Li, Y., and Lence, B. J. (2007a). "Estimating Resilience for Water Resources Systems." Water Resources Research, 43(7), W07422. Li, Y., and Lence, B. J. (2007b). "Mean Points of Failure and Safety Domains: Estimation and Application." Probabilistic Engineering Mechanics, 22(3), 257-266. Li, Y., and Lence, B. J. (2007c). "Numerical Approximations of Design Points in Reliability Analysis under Parametric Changes." ASCE Journal of Engineering Mechanics, in press. Loaiciga, H. A. (2005). "On the Probability of Droughts: The Compound Renewal Model." Water Resources Research, 41(1). Madanat, S. M., and Ibrahim, W. H. W. (1995). "Estimation of Infrastructure Transition Probabilities from Condition Rating Data." ASCE Journal of Infrastructure Systems, 1(2), 120-125. Madsen, H. O., Krenk, S., and Lind, N. C. (1986). Methods of Structural Safety, Prentice-Hall, Englewood Cliffs, NJ. Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (2001). "First-Order Reliability Method for Estimating Reliability, Vulnerability and Resiliency." Water Resources Research, 37(3), 779-790. Melchers, R. E. (1999). Structural Reliability Analysis and Prediction, John Wiley & Sons Ltd., Chichester. Melching, C. S., and Yoon, C. G. (1996). "Key Sources of Uncertainty in QUAL2E Model of Passaic River." ASCE Journal of Water Resources Planning and Management, 122(2), 105-113. Mishalani, R. G., and Madanat, S. M. (2002). "Computation of Infrastructure Transition Probabilities Using Stochastic Duration Models." ASCE Journal of Infrastructure Systems, 8(4), 139-148. North, M. (1980). "Time-Dependent Stochastic Model of Floods." ASCE Journal of the Hydraulics Division, HY5, 649-665. Rice, S. O. (1944). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 23(3), 282-332. Rice, S. O. (1945). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 24(1), 46-156. Salas, J. D. (1996). "Analysis and Modeling of Hydrologic Time Series." Water Resources Handbook, L. W. Mays, ed., McGraw-Hill, 19.01-19.72. 13 Salas, J. D., Chung, C.-H., and Fernandez, B. (2001). "Relating Autocorrelations and Crossing Rates of Continuous- and Discrete-Valued Hydrologic Processes." ASCE Journal of Hydrological Engineering, 6(2), 109-118. Salas, J. D., Delleur, J. W., Yevjevich, V., and Lane, W. L. (1980). Applied Modeling of Hydrological Time Series, Water Resources Publications, Littleton, CO. Sundararajan, C. (1995). "Probabilistic Structural Mechanics Handbook." Chapman & Hall, New York. Tung, Y.-K., and Hathhorn, W. E. (1988). "Probability Distribution for Critical DO Location in Streams." Ecological Modelling, 42(1), 45-60. Tung, Y. K. (1996). "Uncertainty and Reliability Analysis." Water Resources Handbook, L. W. Mays, ed., McGraw-Hill, 7.1-7.65. Tung, Y. K., and Mays, L. W. (1980). "Risk Analysis for Hydraulic Design." ASCE Journal of die Hydraulics Division, 106(HY5), 893-913. Tung, Y. K., Yen, B. C , and Melching, C. S. (2006). Hydrosystems Engineering Reliability Assessment and Risk Analysis, McGraw-Hill, New York. Tyagi, A., and Haan, C. T. (2001). "Uncertainty Analysis Using Corrected First-order Approximation Method." Water Resources Research, 31'(6), 1847-1858. UNESCO/IHP. (1999). Sustainability Criteria for Water Resource Systems, Cambridge University Press, Cambridge, UK. USACE. (2006). "VI-6: Reliability Based Design of Coastal Structures, Coastal Engineering Manual." US Army Corps of Engineers. Weiler, K., Walter, M. T., Walter, M. F., Brooks, E. S., and Scott, C. A. (2000). "Seasonal Risk Analysis for Floodplains in the Delaware River Basin." ASCE Journal of Water Resources Planning and Management, 126(5), 320-329. Yanmaz, A. M. (2000). "Overtopping Risk Assessment in River Diversion Facility Design." Canadian Journal of Civil Engineering, 27(2), 319-326. Yevjevich, V. (1972). Stochastic Processes in Hydrology, Water Resource Publications, Fort Collins, Colorado. 14 C H A P T E R 2 M E A N POINTS O F F A I L U R E A N D S A F E T Y DOMAINS: EST IMAT ION A N D APPLICATION A version of this chapter has been published in the July 2007 issue of Probabilistic Engineering Mechanics. The paper is titled as Mean Points of Failure and Safety Domains: Estimation and Application by Y. Li and B. J. Lence. 15 PREFACE A stochastic process consists of a set of dynamic (or time-dependent) generalizations of variables, each of which is called a sample path or a time function. The system behavior described by a performance function for the system, which may involve several random variables or component processes, is a stochastic process. The nature of two-state transition probabilities is the portion of sample paths of the stochastic process of system behavior that travel from one system state at an earlier time moment into the other system state at a later time moment. Considering this nature, it can be foreseen that knowing the mean points of failure and safety domains may help in estimating two-state lag-1 transition probabilities, and moreover resilience. This chapter defines the mean points, and examines the feasibility and effectiveness of applying them for estimating the transition probabilities. A conceptual representation of the mean points of failure and safety domains in reliability problems is presented, and a practical approach for approximating the mean points based on the first-order reliability method (FORM) is proposed. The effectiveness of the approach is demonstrated with two nonlinear examples. The mean points and the approach for approximating them are then used to develop an innovative method for estimating transition probabilities between failure and safety domains in discrete time-variant reliability problems. The performance of the new method is demonstrated in example problems. It is shown that, from a practical viewpoint, the method is quite effective. Although the example problems are stationary linear and nonlinear univariate systems, the method shows promise for being applicable for more complex cases, such as non-stationary multivariate systems. In the next Chapter, the new method is further generalized to accommodate multivariate cases. 16 2.1 INTRODUCTION In civil engineering reliability analyses, it may be useful to estimate the mean points of the failure and safety domains for a given system. For example, both points may help in estimating two-state transition probabilities between both domains for time-variant reliability problems in the discrete time domain, or their derivatives such as resilience. Estimates of the transition probabilities may be useful in determining the mean number of crossings or the mean interval between crossings in the analysis of dynamic loads on structures, or the likely length of water shortages in water resource allocation strategies. This paper introduces definitions of the mean points of the failure and safety domains, and develops a practical approach for approximating these mean points based on the first-order reliability method (FORM). As one of the potential applications of the mean points, an innovative method of estimating two-state transition probabilities in discrete processes is then proposed based on the mean points. The method is demonstrated for simple cases constructed on a univariate Gaussian (or normal) process. The results of this application indicate that the method is effective in such cases and shows promise for accommodating general cases described by a nonlinear performance function constructed on multivariate and non-normal processes. The remaining sections of this paper are organized as follows. In Section 2.2, the mean points are mathematically defined. In Section 2.3, first-order approximations of the mean points are developed based on FORM, and demonstrated in two nonlinear numerical examples. The application of the mean points for determining two-state transition probabilities is proposed in Section 2.4. Section 2.5 presents conclusions. 2.2 DEFINITIONS Reliability is defined on the basis of a load-resistance interference featured by a performance function g, which involves a set of stochastic variables (x) and a set of deterministic parameters. States of g are divided into two categories, safety (g > 0) or failure (g < 0). The boundary between both states is called the limit-state surface where g =. 0. Thus the probabilistic x-space is divided into two domains by the two states, i.e. the failure domain F where g < 0, and the safety domain S where g > 0. At a given time moment, the probability of a system being in S is referred to as reliability (or instantaneous reliability), while the complement of reliability is known as the failure probability. Practical estimation methods for reliability and failure probability include 1) sampling methods such as those based on Monte Carlo simulations (MCS), for instance Crude Monte Carlo simulations (CMC), 2) approximations such as those based on point estimate, or 17 based on reliability index for instance FORM and the second-order reliability method (SORM), and 3) combinations of above two. (Under the reliability index-based approximations, the reliability index, /?, is the ratio of the mean of g to the standard deviation of g.) For more details regarding reliability methods, see relevant references (Ang and Tang 1984; Madsen et al. 1986; Sundararajan 1995; Tung 1996; Melchers 1999; Lian and Yen 2003). The mean point of failure domain (MPFD), X F , represents the probabilistic geometric centroid of F, mathematically defined as U /(x)dx U /(x)dx x F l = ± = - ^ — (1=1,2,.../.) (2.1) jF/(x)dx PF where x = [ x / . . . x „ ] T , the vector of random variables;/(x) = the joint probability density function (PDF) of x; xpi= the ith coordinate of Xp; and Pf= failure probability. Similarly, the mean point of the safety domain (MPSD), xs, represents the probabilistic geometric centroid of S, mathematically defined as _ j [ * , / ( x ) d x _ £ * f / ( x ) d x £ / ( x ) d j P (i=l,2,...n) (2.2) s where xsi = the ith coordinate of xs; and Ps = reliability. Oftentimes a joint PDF is quite difficult to derive, and the integrals involved in Equations (2.1)-(2.2) may be impossible to solve analytically in complex systems. Therefore, except in rare cases, xs and X F can only be approximated in practice. Numerical schemes and MCS may be used to approximate xs and X F , but their computational costs could be very high. This paper introduces a practical method for approximating xs and X F based on FORM. 2 .3 F I R S T - O R D E R M E A N P O I N T S O F F A I L U R E A N D S A F E T Y D O M A I N S 2.3.1 F O R M Under FORM, the original x-space is transformed into a standardized uncorrelated normal space, u , with zero means and unit standard deviations. To simplify the solution, FORM replaces the original g in the u-space with its first-order Taylor series expanded at the design point, u * as shown in Fig. 2.1. Note that, in this paper, performance functions in both x-space and u-space are identically denoted as g, though they may have different forms. As a local optimum on the limit-state surface of g, u * is locally closest to the origin of the u-space. If g > 0 at the u origin, /? is equal to the distance from the u origin to the limit-state surface, or the norm of u * ; otherwise, 18 f3 is the negative of that distance. Because u * is a local optimum, the following relationship holds: u* = Pa (2.3) where a = the negative of the normalized gradient vector of g at u * . FORM approximates Ps as follows: Ps»&(0) (2.4) where <£( ) = the cumulative distribution function (CDF) of the univariate standard normal distribution. FORM has been shown to be adaptable for complex problems, an advantage over analytical methods, and computationally efficient, an advantage over sampling methods. However, one must be aware of the weakness of FORM. The first-order Taylor series expansion of g may severely distort a highly nonlinear g, reducing its accuracy, and routine optimization algorithms and iterative schemes for finding u * are ineffective for highly irregular g functions, reducing the robustness of FORM for these problems. Hence combinations of FORM and sampling methods are also widely applied in practice in order to improve FORM, particularly importance sampling methods. a first-order Taylor series expansion of g the original g surface safety domain (S) a linearized limit-state surface failure domain (F) a axis Figure 2.1. Schematic description of FORM in a bivariate u-space 2.3.2 First-order approximations of MPSD and M P F D Since the u-space is a standardized uncorrelated normal space, and FORM approximates the limit-state surface as a hyperplane, it is possible to evaluate the analytic solutions to the integrals in Equations (2.1)-(2.2) in the FORM-based u-space, thus producing first-order approximations of the MPFD of the u-space ( u F ) and the MPSD of the u-space ( u s ) , which are denoted as U F * and u s * , respectively. All components of u are standardized uncorrelated normal variables, so 19 U F * and u s * should be located on the a axis. By transforming U F * and u s * back to the x-space, approximations of X F and xs may thereafter be obtained. Because of their central role in approximating X F and x s , u F * and u s * and their approximating approach are discussed herein. 2.3.2 .1 M P F D Since g is approximated as a normal variable in FORM, Equation (2.4) is derived on the basis of a standard normal variable z that is defined as z=*Z*Sl (2.5) o"(g) where E() = the mean, or the expected value; and a() = the standard deviation. Fig. 2.2 shows z and its distribution. Both u * and U F * can be mapped into the z-space, and their analogous values in z-space are denoted as zp and ZF, respectively. In the z-space, ZF can be easily obtained by analytically integrating the numerator in Equation (2.1), as expressed in the following: _ _ Z z ^ z ) d z _ ^ ~ ^ e ^ Z dZ _ -Aizp) _ -<t>A-P) PF Z ^ ( z ) d z * i Z f i ) * { ~ P ) where <fh() = the PDF of the univariate standard normal distribution. Hence in u-space, uF* = -z Fa = ^ < x (2.1) F $>(-P) (a) (b). Figure 2.2. Illustrations of z and its normal distribution (a: in a bivariate u-space; b: in z-space) 20 The result is concise, and straightforward to compute. The properties of up* can be further observed by studying ZF in z-space (Fig. 2.3). It is clear that as zp asymptotically approaches -oo, ZF asymptotically approaches zp, thus UF* can be approximated with u*. However, within the zp range of practical interest, the difference between up* and u* is still considerable. For example, when zp = -3.0,' ZF = -3.2831, the difference being 0.2831. As zp asymptotically approaches +<*>, ZF asymptotically approaches zero, thus U F * can be approximated with the u origin. C c3 r : 1 1 I , - .1 !. .! -•• 1 1 . " ;3 -2.: -1' 0, 1 -<2- 3; Zp Figure 2.3. ZF and zs as function of zp 2.3.2.2 MPSD Likewise, us* can be obtained based on Fig. 2.2, specifically, r^)dz us* = -zsa = =^a (2.8) (2.9) where zs = the analogous value of us* in z-space. As shown in Fig. 2.3, as zp asymptotically approaches -oo, Zs asymptotically approaches zero, thus us* can be approximated with the u 21 origin. As zp asymptotically approaches + ° o , zs asymptotically approaches zp, thus us* can be approximated with u*. 2.3.3 Numerical examples If g is a linear function of u, U F and us simply equal U F * and us*, respectively. For nonlinear gs, the accuracy of u F * and us* varies on a case-by-case basis. Two numerical examples are provided to demonstrate the application of the FORM-based approximation approaches. Here, Equations (2.7) and (2.9) are applied, and CMC is used to generate ten estimates of up and us, each of which samples 2000 realizations. The results of both examples validate the effectiveness of proposed approximations of U F and us. Example 1 The first example is based on the following g function: g =4-0.2244x' + 2 0 X 2 (3-5*' + 1 3 - 5 * 2 r" 6 7 (2.10) x, + x2 4x, + 20x2 where x\ and xi = random variables. Both random variables follow a normal distribution. Their statistical properties are shown in Table 2.1, and the correlation between them is 0.845. The results (Fig. 2.4) indicate that FORM and CMC give similar results for PF, and uF* and u s* are very close to the estimates of U F and us based on CMC. Also, CMC-based estimates of u F experience more scatter than CMC-based estimates of Us. This is an understandable phenomenon, because, in this case, Ps is larger than PF, and as a result, CMC samples fewer realizations in F than in S. Table 2.1. Statistical properties of random variables in numerical examples Example Random variables Mean Standard deviation 1 20 5 x2 16 6 2 Xl . 0.35 /d 0.10/d Xl 0.7/d 0.2/d Xi 16.09 km/d 4.827 km /d x4 18 mg/1 5 mg/1 X5 1.0 mg/1 0.3 mg/1 22 1 — ; 1^ +: uF by CMC • : uF* x: u s by CMC 7 : us* P= 1.1556 PF by FORM = 0.1239 Pr bv CMC = 0.1135 to 0.1265 "-•3.5 0 0.5 1 .. •'. 1.5 .-• 2 Figure 2.4. MPFD and MPSD in numerical example 1 Example 2 The second example examines a more complex g provided by Tung (1990), which is a typical water quality emission response relationship based on the Streeter-Phelps model (Streeter and Phelps 1925) x2 — xl where x5 = mutually independent random variables. Each of these variables follows a lognormal distribution, and their statistical properties are provided in Table 2.1. Fig. 2.5 illustrates five coordinates of up, us, U F * , and us*. The results indicate that FORM slightly overestimates Pp. Most coordinates of the U F * and us* agree closely with those of the u F and us found using CMC, except that the 4th coordinate of U F * is slightly lower than its CMC equivalents. Again, the CMC-based U F estimates experience more scatter than the CMC-based us. 1.5 0.5 s 23 Coordinates of U F , U S , U f * , and Us* Figure 2.5. MPFD and MPSD in numerical example 2 2.4 ESTIMATES OF TWO-STATE TRANSITION PROBABILITIES As one of fundamental properties of reliability problems, MPSD and MPFD may be practically applied in civil engineering reliability analyses. To demonstrate a potential application, a MPSD- or MPFD-based method that estimates transition probabilities between F and S in time-variant reliability problems for a discrete g series is developed herein. (Another potential application of the mean points, which is not further investigated in this paper, could be related to importance sampling. In importance sampling, U F * could be a better candidate than u * for the mean of the sampling density function, because u F * is located inside F, unlike u * that is located on the boundary between F and S. Thus a sampling density function centered on U F * might sample more realizations in F than one centered on u* . ) 2.4.1 Two-state transition probabilities In a discrete g series, the two-state transition probabilities between F and S in two consecutive time steps can be conceptually defined as follows: PFS(t) = P[g(t)<Ong(t+\)>0] (2.12a) PsKt) = P[g(t) > 0 fl g(t+l) < 0] (2.12b) where t = time moment; g(t) = the value of g at t\ g(t+\) - the value of g at r+1; P( ) = probability; PFs(t) = the one-step transition probability from F to 5 at t; and PsiAt) = the one-step transition probability from S to F at t. Note that the time step here represents a time unit defined in a given case, though multiple-step transition probabilities may also be used in practice. (In the 24 continuous time domain, the time step may reduce to an infinitesimal increment.) The transition probabilities in both discrete and continuous time domains hold a central place in time-variant reliability analyses for either non-repairable or repairable systems, and studies regarding these probabilities are also denoted as run or crossing theory. Typical relatives or derivatives of the one-step transition probabilities in both continuous and discrete time domains include various rates or intensities (failure, barrier-, level-, threshold-, scalar-, vector-, up-, down-, or out-crossing), the mean number of crossings, the mean interval between crossings, return period, run (or excursion) length, hazard function, first passage (or excursion, crossing) time, and resilience. Some of these are actually equivalent in nature either unconditionally or conditionally, but are designated differently in different disciplines. One of the cores of these derivatives lies in Equation (2.12) or its analogy in the continuous time domain. For further details on these well-known concepts, see classic descriptions (Rice 1944, 1945; Kapur and Lamberson 1977; Hashimoto et al. 1982; Bras and Rodnguez-Iturbe 1985; Madsen et al. 1986; Melchers 1999). Equation (2.12) represents two intersection probabilities within a stochastic process framework. These probabilities and their derivatives or relatives may be solved through counting synthetic data generated by MCS or a representative historical series (Hashimoto et al. 1982; Tickle and Goulter 1994; Fowler et al. 2003; Karamouz et al. 2003; Luz and Loucks 2003; Kjeldsen and Rosbjerg 2004), but the computational cost of such an approach could be extremely high. In general, analytic solutions to Equation (2.12) only exist in ideally simple systems, e.g. for cases where g is a linear function of a stationary univariate Gaussian process (Sen 1976; Rosbjerg 1977). The analytic solutions are derived on the basis of a bivariate normal distribution (BND) of both g(t) and g(t+\) (Tick and Shaman 1966), and the correlation coefficient between two consecutive gs, pg(\), is readily obtained based on the lag-1 autocorrelation coefficient of the Gaussian process. In complex systems where g is a nonlinear function of multivariate and/or non-normal variables, Equation (2.12) can only be approximated in practice. Primarily there are two practical strategies for approximating /Vs and PSF, one developed based on a BND of both g(t) and g(H-l), and the other on a stochastic linear prediction (SLP) of g(r+l) using g(t). Both have their own shortcomings. For other strategies such as that based on the parallel-system sensitivity measure relative to t, see relevant references (Hagen and Tvedt 1991, 1992). In the BND-based method, all components of x at any time moment are expressed as a linear combination of a given set of standard normal random variables (w) by an appropriate discrete representation approach, such as the Karhunen-Loeve (KL) expansion method (Loeve 1977; 25 Paez and Hunter 2000). The number of components in w depends on the memory of the x process. All gs are thus rebuilt in the w-space by replacing x with w. Equation (2.12) then becomes a parallel-system reliability problem in the w-space. Since Ps can be solved individually by FORM, pg(l) is approximated as a scalar product of the a at t and the a at r+1, which is well-known practice in parallel-system reliability problems (Ang and Tang 1984; Madsen et al. 1986). For further details regarding the BND-based approximation, see relevant references (Der Kiureghian 2000; Maier et al. 2001; Koo and Der Kiureghian 2003; Andrieu-Renaud et al. 2004). However, the implementation of this method may be cumbersome. As previously mentioned, in the representation procedure, the dimension of w depends on the memory of the x process. If an n-dimensional x process has a lag-m covariance structure, w will have (m+\)n components, and after expressing each component of x with w, the rebuilt g will involve (m+\)n random variables in total. For instance, given a lag-1 covariance structure, there are 2n random variables involved in a rebuilt g, as shown in Maier et al. (2001). Remember that only n random variables are involved in the original performance function, and thus the total number of random variables is increased by mn during the representation procedure. As m or n becomes large, the number of random variables in the rebuilt g will increase considerably, making it more difficult to solve Ps at a single moment. In addition, the representation procedure could be complicated. For example, the KL expansion requires eigenvalues and eigenvectors of the autocorrelation matrix of the x process. The SLP-based approximations are mainly typified by those methods that are developed based on the well-known Rice's formula (Rice 1944, 1945). Originally, the formula was derived to estimate the up-crossing rate in a continuous stationary linear univariate Gaussian process. In nature, Rice's formula utilizes a lag-1 SLP of g that requires the time step to be infinitesimal. The SLP is constructed using the derivative process of the original process relative to t. The formula is also approximately applied in a discrete manner. In civil engineering reliability analyses, the Rice's formula-based approximation is mostly applied for estimating PFS and PSF in discrete linear univariate normal cases (Rodriguez-Iturbe 1968, 1969; Nordin and Rosbjerg 1970; Bras and Rodriguez-Iturbe 1985; Salas et al. 2001; Salas et al. 2005), though the BND-based method may provide analytic solutions for these cases. To apply the Rice's formula-based approximation in multivariate cases, the random effects of multiple random variables may be combined into a single variable (Li et al. 2005). However, it has been warned that the discrete application of Rice's formula could cause considerable errors (Nordin and Rosbjerg 1970; Kottegota 1980), or require a complex implementation in practice (Putcha and Patev 1998). 26 2.4.2 A m e t h o d o f e s t i m a t i n g t r a n s i t i o n p r o b a b i l i t i e s b a s e d o n M P S D a n d M P F D MPSD and MPFD provide a new SLP-based method for approximating PFS and PSF- The new method rebuilds the original g, but the rebuilding procedure results in fewer random variables than the existing BND-based method. Developed in the discrete domain, the new method may be more suitable for discrete cases than the Rice's formula-based approximation. As a starting point for the development and application of this new method, its essential strategy is demonstrated here using two stationary univariate examples, one of which is so simple that PFs and PSF can be analytically solved. It is worth stressing that the purpose of testing the examples is merely to demonstrate the essential strategy of the new method, as well as the potential of the new method in more general problems constructed on non-stationary multivariate non-normal processes. 2.4.2.1 A s t a t i o n a r y l i n e a r u n i v a r i a t e e x a m p l e Consider a linear univariate repairable system built on a stationary Gaussian process V where v(t) - the value of V at t; and C = a constant. In addition, assume V to be a first-order autoregressive model, or AR(1), with zero mean and unit variance, as defined in the following: where v(M) = the value of V at r-l; p{\) = the lag-1 autocorrelation coefficient of V; e = a normal white noise with zero mean and unit variance, which is independent of previous vs. See Fig. 2.6 for a schematic of this problem. Since Pps and PSF are similar in concept, the theoretical development here focuses on PFs, without any loss of generality. In this simple example, pg(l) = p(\). At any time moment fo. the analytic solution to Equation (2.12a) is available in this case (Tick and Shaman 1966) g(t) = v(0 - C (2.13) (2.14) P M = £ 2^[v(f0),v(f0+l),/7(l)]dv(f0 + l)dv(r0) (2.15) 27 the PDF of standard Figure 2.6. A stationary Gaussian process V As showed in Fig. 2.6, the probability of v(fo) falling on any point in F at fo, y, can be expressed as <p\{y)dy, where dy is an infinitesimal segment. The probability that v(fo) falls on y and g(fo+l) > 0 equals P[LPx(y) -C> 0]$(y)dy, where LPx(y) = p(\)y + yjl- p(l)2e, a one-step linear prediction of v(fo+l) based on y at to, according to Equation (2.14). Therefore, PFSOO) may be written as follows: PFS CO )= I P[LPx(y) - C>0] </>x (y )dy £p[Z^(y) -C>0]d(y)dy f = EF{P[LPl(y)-C>0]}P[g(t0)<0] (2.16) where Ef( ) = an expectation over F. In order to estimate the expectation, an approximation similar to that employed in developing FORM may be applied, i.e. expanding the probability function in EF as a Taylor series. An important consideration in the expansion is the selection of a proper expansion point. Two candidate points are available, the design point C and the MPFD, Vf(to). The latter is used herein as the expansion point, making the expansion more practical. Using a first-order approximation and noting that Ep(y) = v/r(fo), Equation (2.16) simply reduces to PFS(t()) - P{LPl[vF(t0)]-O0}P[g(t0)<0] (2.17) The physical meaning of Equation (2.17) can be clearly described in terms of sample paths of V (Fig. 2.7). The original state, i.e. Equation (2.12a), include all sample paths from F at fo to S at 28 fo+1, while the approximation by Equation (2.17) considers only sample paths from vF(to) to S at to+l. Here vF(to) is considered to probabilistically represent the entire F domain at fo- It should be noted that the original Rice's formula only takes account of sample paths from point C at fnto S at fn+1 (Melchers 1999). Point C cannot effectively represent the entire F domain, and probably this is the reason why Rice's formula may cause large errors in discrete applications. The approximation based on Equation (2.17) may provide practical convenience in the analysis of PFS, PSF, and their derivatives in complex systems. As an example, the one-step resilience at fo, Re(to), namely a conditional probability that g(fo+l) >0 given g(fo)< 0, can be expressed as P (t ) Re(t0) = F s K o J «F{LP,[v F (f 0 )]-C>0} p[g(f0)<0] | L f V o A I (2.18) A f / * (b) '0+l s l (a) Figure 2.7. Sample paths of V between fo and fo+1 (a: original; b: approximated by Equation (2.17)) In this simple problem, Equations (2.17) and (2.18) can be determined as AA-/>0) R e i t ^ ^ M ^ (2.19) (2.20) Figs. 2.8-2.10 illustrate the results of computational experiments for approximating PFSOO) and Re(to) in the system. It is clear that the approximations by Equation (2.19) properly characterize the overall trend of the analytic solutions given in Equation (2.15) for C < 0. The analytic solution for C = 0 is also available in Salas et al. (2001). When C < -2.0, the difference 29 between the approximation and the exact solution is too small to illustrate schematically. For example, when C = -2.0, the maximum error based on Equation (2.19) is -0.0008, which occurs for p{\) = 0.9. When C < -2.5, Equations (2.19) and (2.15) give the same results as p{\) ranges from -0.9 to 0.9. However, positive C values cause large errors, which are not shown. For comparison, Rice's formula-is also employed in a discrete manner to estimate PFSOO) using the autocovariance function approximated with a finite difference scheme provided by Rodrfguez-Iturbe (1968), as shown in Fig. 2.8. Although the solution accurately converges to the analytic solution as p{\) approaches to 0.9, the discrete application of Rice's formula results in a trend that is biased relative to the analytic solution. As also shown in Nordin and Rosbjerg (1970), the discrete application of Rice's formula poorly estimates /Vs(£o) when both C and p{\) are zero, which is a quite typical scenario examined in stochastic process modeling. 0.5, r —1 1 I i i i L_ -0.9 -0.6 -0.3 0 0.3 0.6 0.9 Figure 2.8. Computational experiments on PFsik) (C - 0) 1 * ~ " " - * - - ~ ~ ! — — - ^ ~ - » ^ \AC = -0.5 C = 0I*" f~--^_ \ ^**->>^ y N ^ x V i The analytic solution —"—The approximation based on Equation (2.20) | -i i i i i i 0 I i 1 1 i i i 1— -0.9 -0.6 -0.3 0 0.3 0.6 0.9 AD Figure 2.10. Computational experiments on Re(to) (C = 0 to -3.0) This simple example is a widely studied topic in engineering practice. At any time moment, the example actually represents a special interpretation of FORM in that u * = C and ft - -C. Observations obtained from the computational experiments will be informative for practical applications of the MPSD- or MPFD-based method in estimating two-state transition probabilities. Three observations seem rather important in a practical sense. As mentioned, the maximum error of the proposed method generally increases as C increases, particularly when C > 0. This phenomenon implies that the second-order and higher-order terms in the Taylor expansion series of Equation (2.16) become more significant when C becomes large, for example as C > 0. (Conversely, it may be expected that a negative C will cause a large error in applying the MPSD-based method for estimating PSF-) According to FORM, however, a larger C (or a smaller fj) in this example indicates a more risky system. In reality, the reliability of an important system is normally kept high; e.g., a Ps of 0.9 requires that f3 > 1.3, and a Ps of 0.8 requires that /3 > 0.85. Thus a positive C is basically meaningless from a practical viewpoint. Even if a positive C occurs, the MPFD-based method for approximating Pfs(fo) and Re(to) is still applicable, considering an equivalence of PFsik) resulting from Equation (2.15) Pn Co) = £ £" & M'o X v(/0 + l),p(\)]dv(t0 + l)dv(f0) = JT C & W M f o + D,/?(l)]dv(f0 +l)dv(r0) (2.21) Changing the sign of C to be negative and then re-estimating PFs(to) with -C does not change the final result. Provided that -C is not too close to zero, the MPFD-based method in general gives good results, as shown in Figs. 2.8-2.10. 31 It is also shown that errors increase as the absolute value of /XI) becomes large, for example as p(l) approaches 0.9 in most cases. Fortunately, high absolute values of p{\) close to 1.0 are very unusual in stationary time series for many real engineering applications. In hydrological practice, p{\) is generally small for natural stationary series, annual or monthly, though high absolute values of p{\) may result from special effects such as non-stationarities (Salas 1996). In particular, the noise series extracted by time series modeling oftentimes have moderate values of p{\). Hence the practical applicability of MPSD- or MPFD-based methods is not substantially narrowed. The MPFD-based method demonstrates excellent performance for negative p{\), except for the case with C - 0. This implies that the proposed method might be especially suitable for negatively correlative processes, as long as C is not too close to zero. In general, the MPFD-based method demonstrates satisfactory modeling results for P F S and Re in a fairly typical example, from a practical viewpoint. Considering that the example is a special interpretation of FORM, the results of numerical experiments for the example preliminarily indicate the applicability of MPSD- or MPFD-based methods for approximating PFS, PSF, and Re in stationary, nonlinear, multivariate and non-Gaussian systems. Adverse scenarios might potentially occur when pg(\) is too close to the extreme values, for example near 0.9 or -0.9. 2.4.2.2 A stationary nonlinear univariate example As a special interpretation of FORM, the previous example preliminarily shows the adaptability of the MPFD-based method in more complex systems. Now let's extend the MPFD-based method to a slightly more complex nonlinear example, g(t) = v(t) + em+C (2.22) where v(t) and C have been defined in 4.2.1. This example mainly investigates the ability of the MPFD-based method to solve PFs-It has been shown in the previous example that, the higher the system reliability, the more accurate the MPFD-based method will be, and that extreme lag-1 correlation coefficients of g may cause large errors in the method. To examine the MPFD-based method reasonably, C will be adjusted to maintain moderate system reliabilities, and moderate positive values of p{\) will be used. Usually moderate values of p\\) and PS are typical of normal civil engineering systems. 32 In this nonlinear example, the analytic solution may be difficult to derive, so a CMC simulation with 10000 realizations is employed to produce the "exact" solutions of Psand PFs. The MPFD-based approximation is implemented in the following procedure. First, PF(to) is estimated with the original g using FORM, and uF* is estimated by Equation (2.7). In this example, uF* equals the first-order approximation of the MPFD of the v-space, vF*. Then, according to Equation (2.14), a SLP of v(t0+l) based on vF* is obtained as LPl(vF*) = p(l)vF* + yl\-p(l)2e (2.23) Thus g(?o+l) is rebuilt by replacing v(t0+l) with LP\{vF*), g(t0 +1) = p(\)vF * + P(V2£ + ep0)v"'+^-^E + C (2.24) A reliability estimate with the rebuilt g(to+\) is solved using FORM. Note that this reliability estimate is different from the reliability of this system at ro+1. Finally, PFs(to) is obtained as a product of this reliability estimate and PF(to), which is similar to Equation (2.17) PM = P[pi\)vP * + >/l - Pi}f £ + epm^+^£ +C> 0]PF(tQ) (2.25) It can be expected that, given an n-dimensional x process with a lag-m covariance structure, the rebuilt g of the MPFD-method will involve mn random variables, and the rebuilding procedure increases the number of random variables by (m-l)n. Specifically, in this example, Equation (2.24) involves one random variable, e, and no increase in the number of random variables occurs. For comparison purposes, the BND-based approximation in conjunction with the KL expansion is applied. Since V is an AR(1) process, w in the KL expansion will have two components. Hence the rebuilt g involves two random variables. For C = 2.0, 1.0, and 0, the FORM estimates for Ps are 0.9830, 0.8995, and 0.7174, respectively, while the CMC estimates for Ps are 0.9823, 0.8998, and 0.7186, respectively. Table 2.2 summarizes the computational results for PFs(to) under a range of moderate values of p(\). In this stationary nonlinear univariate example, FORM effectively estimates Ps. The MPFD-based method and the BND-based method produce very close estimates for PFs, which are quite satisfactory results compared with CMC estimates. Considering that the BND-based approximation involves more random variables, however, the MPFD-based method appears simpler to implement in principle. Again, it is found that, either a higher p{\) or a lower Ps will generally cause a larger error in the MPFD-based method. This finding is consistent with that observed in the previous example. 33 Table 2.2. Computational results of PFS in the nonlinear example c Methods AD 0.6 0.7 0.8 MPFD 0.0134 0.0120 0.0100 2.0 BND 0.0132 0.0118 0.0100 CMC 0.0137 0.0132 0.0100 MPFD 0.0615 0.0532 0.0421 1.0 BND 0.0612 0.0534 0.0438 CMC 0.0618 0.0540 0.0429 MPFD 0.1218 0.1012 0.0743 0 BND 0.1241 0.1068 0.0867 CMC 0.1234 0.1034 0.0813 2.4.3 Summary In order to approximate two-state transition probabilities in complex discrete dynamic reliability problems, the performance function may have to be rebuilt within a stochastic process framework. The rebuilding process may change the number of random variables in the resulting performance function. The MPSD- and MPFD-based approximating method leads to fewer random variables than the current BND-based method. In general, a performance function with fewer random variables will be easier to solve in practice. Therefore the MPSD- and MPFD-based approximation could be a promising method in discrete dynamic reliability problems. Further studies are needed to test the performance of the MPSD- and MPFD-based approximation in non-stationary scenarios. Moreover, with a linear forecasting model for all random variables involved, it may become possible to estimate multiple-step two-state transition probabilities using the MPSD- or MPFD-based methods. One should also realize that the feasibility of the MPSD- or MPFD-based methods in a system depends on two important conditions; that is, 1) FORM is sufficiently effective for the system, and 2) a linear forecasting model is available. Whether the former is met is case-specific. The latter is often easily met, because any time-variant reliability problem requires a well-built stochastic process model for all random variables as a prerequisite, such as a time series model; otherwise the problem cannot be properly analyzed. 34 2.5 CONCLUSIONS The MPSD and MPFD have been defined. Both points belong to the fundamental properties of reliability analysis. A FORM-based approximation method for both points has been developed. The effectiveness of the method has been verified using two nonlinear examples. The current BND- and SLP-based approximations for estimating two-state transition probabilities between F and S in discrete time-variant reliability problems may be deficient. This paper develops an alternative method based on MPSD and MPFD for estimating the transition probabilities in discrete processes. The method presents a new way of studying dynamic reliability problems. The central idea of the proposed method is demonstrated in univariate linear and nonlinear examples, and it provides satisfactory results in a practical sense. The new method requires fewer random variables in rebuilding the original performance function than the current BND-based method. Moreover, the new method is especially suitable for discrete processes, unlike the Rice's formula-based method. It is expected that the method can be extended to cases with non-stationary, nonlinear, multivariate and non-Gaussian processes, and this extension, namely a generalized form of the new method, will be investigated in the future work. In addition, other potential applications of MPSD and MPFD, such as for improving importance sampling, will be investigated. 35 References Andrieu-Renaud, C , Sudret, B., and Lemaire, M. (2004). "The PHI2 method: a way to compute time-variant reliability." Reliability Engineering & System Safety, 84(1), 75-86. Ang, A. H. S., and Tang, W. H. (1984). Probability Concepts in Engineering Planning and Design, Vol. II: Decision, Risk, and Reliability, John Wiley, New York. Bras, R. L., and Rodrfguez-Iturbe, I. (1985). Random Functions and Hydrology, Addison-Wesley Publishing Company, Reading, MA. Der Kiureghian, A. (2000). "The Geometry of Random Vibrations and Solutions by FORM and SORM." Probabilistic Engineering Mechanics, 15(1), 81-90. Fowler, H. J., Kilsby, C. G., and O'Connell, P. E. (2003). "Modeling the Impacts of Climate Change and Variability on The Reliability, Resilience, and Vulnerability of A Water Resource System." Water Resources Research, 39(8), 1222. Hagen, 0., and Tvedt, L. (1991). "Vector Process Out-Crossing as Parallel System Sensitivity Measure." ASCE Journal of Engineering Mechanics, 117(10), 2201-2220. Hagen, 0., and Tvedt, L. (1992). "Parallel System Approach for Vector Out-Crossing." ASME Journal of Offshore Mechanics and Arctic Engineering 114(2), 122-128. 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. Kapur, K. C , and Lamberson, L. R. (1977). Reliability in Engineering Design, John Wiley & Sons, Inc. Karamouz, M., Szidarovszky, F., and Zahraie, B. (2003). Water Resources Systems Analysis, CRC Press LLC, Boca Raton, FL. Kjeldsen, T. R., and Rosbjerg, D. (2004). "Choice of Reliability, Resilience and Vulnerability Estimators for Risk Assessments of Water Resources Systems." Hydrological Sciences Journal, 49(5), 755-767. Koo, H., and Der Kiureghian, A. (2003). "FORM, SORM, and Simulation Techniques for Nonlinear Random Vibrations." UCB/SEMM-2003/01, Department of Civil and Environmental Engineering, University of California, Berkeley, California. Kottegota, N. T. (1980). Stochastic Water Resources Technology, The Macmillan Press Ltd., Hong Kong. Li, C. Q., Lawanwisut, W., and Zheng, J. J. (2005). "Time-Dependent Reliability Method to Assess the Serviceability of Corrosion-Affected Concrete Structures." ASCE Journal of Structural Engineering, 131(11), 1674-1680. Lian, Y., and Yen, B. C. (2003). "Comparison of Risk Calculation Methods for a Culvert." ASCE Journal of Hydraulic Engineering, 129(2), 140-152. Loeve, M. (1977). Probability Theory I, Springer, New York, NY. Luz, L. D., and Loucks, D. P. (2003). "Developing Habitat Suitability Criteria for Water Management: A Case Study." International Journal of River Basin Management, 1(4), 283-295. Madsen, H. O., Krenk, S., and Lind, N. C. (1986). Methods of Structural Safety, Prentice-Hall, Englewood Cliffs, NJ. Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (2001). "First-Order Reliability Method for Estimating Reliability, Vulnerability and Resiliency." Water Resources Research, 37(3), 779-790. Melchers, R. E. (1999). Structural Reliability Analysis and Prediction, John Wiley & Sons Ltd., Chichester. 36 Nordin, C. F., and Rosbjerg, D. (1970). "Applications of Crossing Theory in Hydrology." Bulletin of the International Association of Scientific Hydrology, 15(1), 27-43. Paez, T. L., and Hunter, N. F. (2000). "Representation of Random Shock via the Karhunen Loeve Expansion." the 2000 Shock and Vibration Symposium, Washington, DC. Putcha, C. S., and Patev, R. C. (1998). "Time-Variant Reliability Assessment." Uncertainty Modeling and Analysis in Civil Engineering, B. M. Ayyub, ed., CRC Press LLC, Boca Raton, FL, 73-97. Rice, S. O. (1944). "Mathematical Ana lys is of Random Noise." Bell System Technical Journal, 23(3), 282-332. Rice, S. O. (1945). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 24(1), 46-156. Rodriguez-Iturbe, I. (1968). "A Modern Statistical Study of Monthly Levels of The Orinoco River." Bulletin of the International Association of Scientific Hydrology, 13(4), 25-40. Rodriguez-Iturbe, I. (1969). "Applications of The Theory of Runs to Hydrology." Water Resources Research, 5(6), 1422-1426. Rosbjerg, D. (1977). "Crossing and Extremes in Dependent Annual Series." Nordic Hydrology, 8(5), 257-266. Salas, J. D. (1996). "Analysis and Modeling of Hydrologic Time Series." Water Resources Handbook, L. W. Mays, ed., McGraw-Hill, 19.01-19.72; Salas, J. D., Chung, C.-H., and Cancelliere, A. (2005). "Correlations and Crossing Rates of PeriodicTStochastic Hydrologic Processes." ASCE Journal of Hydrological Engineering, 10(4), 278-287. Salas, J. D., Chung, C.-H., and Fernandez, B. (2001). "Relating Autocorrelations and Crossing Rates of Continuous- and Discrete-Valued Hydrologic Processes." ASCE Journal of Hydrological Engineering, 6(2), 109-118. Sen, Z. (1976). "Wet and Dry Periods of Annual Flow Series." ASCE Journal of the Hydraulics Division, 102(HY10), 1503-1514. Streeter, H. W., and Phelps, E. B. (1925). "A Study of The Pollution and Purification of The Ohio River, III, Factors Concerned in The Phenomena of Oxidation and Reaeration." Public Health Bulletin 146, U.S. Public Health Service, Washington, D.C. Sundararajan, C. (1995). "Probabilistic Structural Mechanics Handbook." Chapman & Hall, New York . Tick, L. J., and Shaman, P. (1966). "Sampling Rates and Appearance of Stationary Gaussian Processes." Technometrics, 8(1), 91-106. Tickle, K. S., and Goulter, I. C. (1994). "Statistical Properties of Reliability and Resilience Measures." Stochastic and Statistical Methods in Hydrology and Environmental Engineering, Volume 4: Effective Environmental Management for Sustanable Development, K. W. Hipel and L. P. Fang, eds., Kluwer Academic Publishers, 209-220. Tung, Y . -K. (1990). "Evaluating the Probability of Violating Dissolved Oxygen Standard." Ecological Modelling, 51(3-4), 193-204. Tung, Y . K. (1996). "Uncertainty and Reliability Analysis." Water Resources Handbook, L. W. Mays, ed., McGraw-Hill, 7.1-7.65. 37 C H A P T E R 3 EST IMATING RESIL IENCE F O R W A T E R R E S O U R C E S S Y S T E M S A version of this chapter has been published in Water Resources Research. The paper is titled Estimating Resilience for Water Resources Systems by Y. Li and B. J. Lence. P R E F A C E Chapter 2 defines the mean points of the failure and safety domains, which are identified as a new fundamental property of reliability problems. A stochastic linear prediction (SLP) of the performance function based on the mean points leads to an innovative methodology for estimating lag-1 two-state transition probabilities in discrete processes, and furthermore lag-1 resilience. From a practical point of view, the methodology shows satisfactory modeling effects in systems constructed on stationary univariate normal processes. Identifying whether the methodology is suitable for systems constructed on non-stationary maltivariate non-normal processes is one focus of this chapter. A generalized form of the mean points-based methodology for estimating the two-state transition probabilities and resultant resilience in the multivariate scenario is formulated. This chapter also develops a practical approximation for the autocorrelation coefficient of performance function process, which leads to a new bivariate normal distribution (BND) based methodology for estimating the two-state transition probabilities. Both methodologies are implemented in a water supply example. 39 3.1 INTRODUCTION In civil engineering risk analyses, water resources systems are generally classified into two categories, repairable or non-repairable. Non-repairable water resources systems are typified by those hydraulic facilities whose major failure mode is serious structural failure, such as a dam breach or a levee slide. Once such a structural failure has occurred, these facilities cannot be easily restored to a functional state, and in a practical sense, the failure is considered permanent. In repairable systems, the major failure mode is performance failure, such as in the case of water quality violations, drought or flood events, or water supply shortages. A failed repairable water resources system may return to a functional state, as the system's natural input conditions, hydrological or meteorological, follow the natural cycle. Obviously, the recovery capacity from a system failure to a functional state is an important consideration in the planning, design, and maintenance of repairable water resources systems. In general, this capacity can be measured by resilience. A concept initially developed in ecological science, resilience indicates the ability of a system to persist in the face of forcing factors (Holling 1973; Fiering 1976). By resilience, an ecological system may absorb certain forcing factors and survive, and then continue to function. The concept has also been operationalized in water resources engineering as a measure of the system's recovery capacity. Despite many interpretations of this ecological concept in the water resources field (Fiering 1982), resilience is generally defined in practice as the conditional probability that, given a system failure at a specific moment, the system will recover in a following moment. This definition has been widely accepted in measuring the sustainability of a water resources system (Loucks 1997; ASCE Task Committee on Sustainability Criteria 1998; UNESCO/IHP 1999). Since most systems are designed to function for a given life period, the conditional probability must be estimated within the framework of a stochastic process. Current estimation methods are either computationally expensive or theoretically deficient, as is addressed in Section 3.2. Thus it is important to explore efficient and reliable methods for estimating resilience. This paper proposes two practical methods for estimating the lag-1 resilience in the discrete time domain that are applicable under either stationary or non-stationary conditions. One method is developed based on a bivariate normal distribution (BND) of the performance functions at two consecutive moments, while the other is based on a stochastic linear prediction (SLP) of the performance function using the mean point of the failure domain (MPFD). An approximation method for estimating the MPFD has been developed by Li and Lence (2007). 40 The remainder of this paper is organized as follows. Section 3.2 reviews background information regarding resilience. In Section 3.3, the two aforementioned methods are developed in detail. In Section 3.4, these methods are illustrated with a hypothetical water supply problem of providing water to a community from two sources, in which the demand flows and available supply flows vary stochastically. Section 3.5 presents the conclusions. 3.2 B A C K G R O U N D 3.2.1 D e f i n i t i o n s The foundation of probabilistic indicators such as resilience is a load-resistance interference depicted by a performance function, g, of random variables and deterministic parameters. The probabilistic space of the random variables is divided into two domains, i.e., the failure domain, F, where g < 0, and the safe (or success) domain, S, where g > 0. The boundary between both domains is the limit-state surface where g = 0. The probability of the random variables of a system falling into S is called reliability, while the complement of reliability is the failure probability. Practical estimation approaches for reliability include 1) sampling methods such as Monte Carlo Simulation (MCS), 2) approximations such as point estimation methods, as well as reliability index-based methods typified by the first-order reliability method (FORM), and 3) combinations of the previous two such as importance sampling methods. For further details on reliability methods, please see Ang and Tang (1984), Madsen et al. (1986), Harr (1987), Sundararajan (1995), Melchers (1999), Haan (2002), and Tung et al. (2006). Considering its importance in this paper, FORM is briefly described here. FORM first transforms the original random variables (x) into a standardized uncorrelated normal space (u ) with unit standard deviations and zero means using their statistical distributions, and then replaces the original g with its first-order Taylor series expanded at the design point ( u * ) . A local optimum on the limit-state surface, u * is locally closest to the origin of the u-space. If g > 0 at the u origin, the reliability index (ft) equals the distance from the origin to u * ; otherwise, /? is the negative of that distance. FORM approximates reliability as PS~®(P) (3.1) where Ps = reliability; and $( ) = the cumulative distribution function (CDF) of the standard normal distribution. FORM has been shown in many hydrotechnical applications to be adaptable to complex problems and to be computationally efficient. Nevertheless, its accuracy may be poor in cases where g is highly non-linear, and problems with many random variables might reduce its efficiency. 41 Given a time moment t\ and a later moment t2, resilience is defined as Re(tl,t2) = P[g(t2)>0/g(t])<0] _P[g(tl)<Ong(t2)>0] P[g(tl)<0] (3.2) where Re(t\, t2) = the resilience between t\ and t2, or the conditional probability that, given a system failure at t\, the system will recover at t2; g(t\) and g(t2) = the performance functions at fi and t2, respectively; P[ ] = probability; PF() = the failure probability - and PFs(ti, h) = a two-state transition probability, or the intersection probability that a system fails at t\, and recovers at t2. One should note that the time domain in Equation (3.2) may be either discrete or continuous, and the time lag At (= t2- t\) may be either infinitesimal or considerable. In the discrete time domain, At is replaced with a given number of time steps. The most typical lag is one step, thus Re(t\, t2) reduces to the so-called lag-1 resilience at./i, Re\(t\), which has significant practical applications in stochastic hydrology. If the g of a system is a stationary process, Re\(t\) becomes the reciprocal of the mean duration of a failure event, which is by definition simply the return period of a success event (Rodriguez-Iturbe 1968, 1969; Sen 1976; Rosbjerg 1977; Hashimoto et al. 1982). Moreover, when g is time-independent, the system reduces to a binominal process, also known as a Bernoulli trial; that is, g is an independent and identically distributed (iid) variable. The binominal process forms the essential foundation for common frequency analyses; e.g., an annual maximum stream flow sequence may be typically assumed to be iid. 3.2.2 Previous work The history of developing approaches for evaluating PFSQI, h) and its derivatives including resilience dates back to the 1940s. Initially, the instantaneous crossing rate of a continuous univariate stationary Gaussian (or normal) process over a specific threshold level was studied, leading to Rice's classic formula (Kac 1943; Rice 1944, 1945). At present, three methodologies are widely used in practice for estimating PFs(h, t2) or its derivatives. Summaries of various methodologies for estimating PFs(h, h) or its derivatives are provided by Andrieu-Renaud et al. (2004) and Li and Lence (2007). The first one adopts sampling methods by counting failures and successes among a set of synthetic data series (or ensembles) generated by a time series model for the random variables using MCS. For example, Re\(t\) is estimated in this way as follows: 42 Reit)„Ml£L (3.3) where N/(ti) = the number of failure realizations in the g(ti) ensemble; and N/S(t\, ti) = the number of success realizations in the g(ti) ensemble which evolve from the failure realizations in the g(t\) ensemble. Variance reduction techniques, such as importance sampling, might improve the efficiency of this estimation procedure, but a general strategy for such an improvement is not available at the moment. Under stationary conditions, the g series is usually assumed to be ergodic, thus Re\(t\) can be estimated using one g series that may be a synthetic sequence generated by M C S , or simply a representative historical record (Hashimoto et al. 1982; Fowler et al. 2003; Karamouz et al. 2003; Luz and Loucks 2003; Kjeldsen and Rosbjerg 2004). Equation (3.3) may require an extraordinarily large ensemble to achieve high accuracy. Since PF of a system in normal operations is small, Nj(t\) is relatively small. On the other hand, in order to guarantee an adequately large sample size in evaluating Equation (3.3), a large Nj(ti) is preferable, requiring in turn extraordinarily large ensembles. Assuming Rei{t\) to be 0.5, for example, Nj(t\) should preferably be around 386 to limit the relative error of the 95% confidence interval for estimates of Equation (3.3) to be 10%, based on Equation (3.27) in Section 3.4. Letting F/r(fi) be a quite moderate value, say 0.05, an ensemble size of roughly 7720 wil l , thus, be needed to obtain a Nj{t\) of 386. For M C S estimates of this Pp(ti), an ensemble size of about 7300 limits the relative error of the 95% confidence interval to be 10%, according to Shooman (1968). In this case, therefore, resilience estimation requires an even larger ensemble size than reliability estimation does. Thus the sampling methods might be computationally expensive. The second methodology adopts analytic methods based on a SLP of g, which are typically represented by Rice's formula-based methods. The formula applies to cases where g is a linear function of a continuous univariate stationary Gaussian process. The nature of the formula is a S L P of g{ti) using g(t\) that requires At to be infinitesimal. In stochastic hydrology, the formula has been applied for modeling univariate hydrological time series, such as a water level or stream flow rate (Rodriguez-Iturbe 1968, 1969; Nordin and Rosbjerg 1970; Bras and Rodriguez-Iturbe 1985; Salas et al. 2001; Salas et al. 2005). Although actually continuous in time, most hydrological processes in practice are described in a discrete time scale. Hydrological data are recorded as discrete observations; that is, rain gauges may record rainfall intensities in five-minute intervals. Also, these data are used in a discrete format; e.g., various time-averaged values (daily, weekly, monthly, or annually) are widely used in describing stream flow rates. In the aforementioned hydrological applications, Rice's formula was simplified in a discrete 43 manner. Due to this simpl i f icat ion, however, the formula may provide poor results in stochastic hydrological model ing, and, thus, may be an inappropriate tool for describing hydrological processes (Nordin and Rosbjerg 1970; Kottegota 1980). L i and Lence (2006) examined the abil i ty of R i ce ' s formula to estimate the up-crossing probabil i ty in discrete stationary univariate Gaussian processes, and found that R i ce ' s formula performs poorly in important scenarios, especial ly for estimating resil ience. Addi t iona l ly , R i ce ' s formula may be hard to implement in practice (Putcha and Patev 1998; L i and Lence 2006). See L i et al. (2005) for the complicated analytic solution of the formula. The third methodology adopts analytic methods based on a B N D of g(ti) and g(t\). The key to this idea is the evaluation of the correlation coefficient between g{ti) and g(t\). This methodology results in analytical solutions to PFs(t\, ti) for discrete stationary linear univariate Gaussian processes (Tick and Shaman 1966; Sen 1976; Rosbjerg 1977; Salas et al. 2001; Salas et al. 2005). In complex scenarios with a non-linear g of multivariate stationary random variables, the methodology may be implemented in an approximate manner. Current B N D techniques approximate the correlation coefficient by applying F O R M at both h and t\ (Der K iureghian 2000; Ma ie r et al. 2001; K o o and Der K iureghian 2003; Andr ieu-Renaud et al . 2004). However , these techniques may signif icantly increase the number of random variables, thus compl icat ing the solution procedure for the correlation coeff icient (L i and Lence 2007). T o date, resil ience has been estimated pr imari ly under stationary condit ions; that is, all random inputs to a system are actually or assumedly stationary. Due to natural and art i f icial effects, however, hydrometric data in many cases exhibit various non-stationarities, such as seasonal periodicit ies and annual trends. These non-stationarities certainly inf luence short-term operations and long-term planning (Deo et al. 1997; We i le r et al . 2000; Durrans et al . 2003; Coles and Tawn 2005). One adaptation to non-stationarities divides the time hor izon roughly into several separate segments within each of which all random variables are assumed to be stationary, and then evaluates resil ience within each segment under stationary condit ions (Fowler et al. 2003). However , this approach is unable to represent the continuous temporal variation o f resil ience, wh ich is quite crucial for short-term seasonal operations in many systems. 3.3 PROPOSED METHODOLOGIES This paper proposes two methods for approximating lag-1 resil ience under either stationary or non-stationary condit ions, both of wh ich are applicable to multivariate non-linear, non-Guassian g, without loss of generality. In both methods, a stochastic process is assumed to be 44 able to effectively describe the inner driving mechanism of all random variables. Both methods have two fundamental bases. Firstly, a periodic vector (or multivariate) autoregressive moving-average (VARMA) model, which is a kind of time series model, mimics the stochastic process. For details regarding periodic VARMA models and their applications in stochastic hydrology, see Salas et al. (1980), Kottegota (1980), Salas (1996), Brockwell and Davis (1996), and Chiang and Gates (2004). Secondly, FORM estimates all probabilistic contents required for estimating lag-1 resilience. Major considerations for choosing FORM are its relatively high computational efficiency compared with sampling methods, and its high accuracy compared with other simpler reliability index-based approximations, such as the mean value first-order second-moment (MVFOSM) method. The basic information regarding periodic VARMA models used in developing both methods is herein briefly described. A VARMA model decomposes the data series of all random variables involved in a system into deterministic components and stationary random residuals. The data series may be original, or transformed if necessary. The deterministic components include periodic (or seasonal) and long-term trends. The random residuals, also ergodic, form a set of random variables to be investigated with probabilistic analyses, x = (xi, x2,..., x„)T, which is transformed to a set of standardized random variables, u = u2,..., w„)T. Here T denotes the transpose. Both x and u are time-dependent; that is, x = x(f) and u = u(f) where t is time. In addition, the first-order dependence structures of x are shown as follows: 1 xO rx l = M U ) V^o(">l) 7xM,n) (3.4a) (3.4b) where T\Q = the lag-0 covariance matrix of x; I"xi = the lag-1 covariance matrix of x; %o(i,j) = the covariance between x,(f) and x/r); and yx\(i,j) = the covariance between x,(f + 1) and x/r)- It should be noted that building a periodic VARMA model, particularly covering the whole year, usually demands intensive supporting data. 3.3.1 T h e b i v a r i a t e n o r m a l d i s t r i b u t i o n ( B N D ) m e t h o d In FORM, g is approximated as a normal random variable, and, thus, is represented by a standardized normal variable z 45 Z = 1 ^ M 3 5 ) °(g) where E() = mean, or expected value; and c() = standard deviation. In a discrete time domain, values of z at two consecutive moments in time, t\ and ti (= t\ + At), follow a BND. Clearly, PFs(t\, h) is the probabilistic content of the shaded area in Figure 3.1, as defined in the following: PFS - h) = C' {[z(f|}' z(t>X p[z(ti}'z(f2)] ]dz(t>)ck(f>} (3 '6) where z(t\) and z(?2) = the standardized normal variables for g(t\) and g{ti), respectively; p\t\) and p\t2) = the reliability indices at t{ and r2, respectively; p[z(t\), z{h)] = the correlation coefficient between z{t\) and z(r2); and 0j() = the probability density function (PDF) of the BND. ////, Z(t2) • -m 0 z(ti) - /*(«) Figure 3.1. A probabilistic space of z(t\) and z(?2) By definition, p[z(t\), z(t2)\ is expressed as Cov[z(f , ) ,z(f 2)] Piz(tM(h)>~:TTr7* ( 3 J ) o[z(r,)]a[z(r2)] where Cov( ) - covariance. Both z(t\) and z{h) are standard normal variables, so their standard deviations are one. Substituting Equation (3.5) into Equation (3.7), therefore, ^ , ) , z ( r 2 ) ] = C ° v [ g ( f | ) ' g ( f ' ) ] (3.8) o-[gO-,)Mg(f2)] To approximate the denominator of Equation (3.8), g can be replaced with its first-order Taylor series expanded at u* g » g(u*) + Ju T • (u - u*) (3.9) where Ju = the gradient vector (or Jacobian) of g relative to u at u*. Noting that u is a standardized normal space, 46 « t ) A & f H % - ) ' - W (3.10) ^ ow, au2 aun Thus a[g(0]-||Ju(0| (3.11a) o[g(r2)] = |ju(r2)| (3.11b) where Ju(r() and Ju(r2) = the Jacobians relative to u at the design point for t\ and t2, respectively. Approximating the numerator of Equation (3.8) requires knowledge of the lag-1 dependence structure of u. Although Txi is typically estimated based on statistical analyses of original data series, it is difficult to develop the dependence structure of u in reality due to the transformation from x to u. To overcome this difficulty, g is expanded at x*, the mapping of u* in x-space. Ignoring second-order or higher-order terms, g reduces to g = g(x*) + JxT-(x-x*) (3.12) where Jx = the Jacobian of g relative to x at x*. By manipulating the covariance on two such expansions at t\ and t2, one obtains C o v ^ O . ^ l - J x ^ - r x i - J x a , ) (3.13) where Jx(fi) and Jx(f2) = the Jacobians relative to x at x* for t\ and t2, respectively. Hence Equation (3.7) is expressed as a i 4 ) Knowing p[z(t\), z(t2)], p\t\), and p\t2), PFSOI, h) can be obtained using Equation (3.6). Finally, Re\(t\) is estimated using Equation (3.2) in which /V(?i) can be estimated with any appropriate approach such as FORM or sampling methods. A special scenario is noted. When t\ = t2, x(fi) and x(f2) come from the same probabilistic space; that is, they are exactly same. Then p[z(ti), z(t2)] becomes the correlation coefficient between two different g in one probabilistic space. In this case, Equation (3.9) can be conveniently used for estimating the covariance in Equation (3.8), considering the lag-0 covariance matrix of u in fact is a unit matrix. For example, assuming two performance functions atf,,gi(fi)and g2(tx), Cov[s,(r,), g2(tl)l = Ju,^) 1 • Ju2(/,) (3.15) where Jui(fi) and Ju2(fi) = the Jacobians of gi(fi) and g2(fi) relative to u at their design points, respectively. Then Equation (3.8) becomes 47 p U ( t ) z ( t ) 1 = J u , ( ? , ) T - J u 2 ( r , ) Piz^z^)\ | J U ] ( r i ) | | | j U 2 ( f | ) | = a,(0T-a2(0 (3.16) where Z\(t\) and z2(t\) = the standardized normal variables for g\(t\) and g2(ti), respectively; and rxi(fi) and Gl2(?i) = the negatives of the normal ized gradient vector of gi(fi) and § 2(? 1) at their design points, respectively. Equat ion (3.16) is s imply the wel l -known F O R M - b a s e d approximation for the correlation coefficient between two failure modes in system rel iabi l i ty problems (Ang and Tang 1984; Madsen et al. 1986). Equat ion (3.16) is only applicable when all random variables come f rom one probabil ist ic space. The misuse of Equat ion (3.16) for estimating Re\(t[) and Pps(ti, t2) may cause misleading results, because, according to the principles of stochastic processes, x(ti) and x(t2), values of a process at different time steps are different random variables f rom two different probabil ist ic spaces. The fo l lowing example demonstrates the potential results of misusing Equat ion (3.16). Assume that g(t[) and g(t2) have the same form with \(ti) and x(t2) f rom two different spaces. Under stationary condit ions, g(t[) and g(t2) have the same values for their design points, rel iabi l i ty indexes, and normal ized gradient vectors. Thus Equat ion (3.16) results in p[z(t\), z(t2)] = 1, and Tx i is, thus, neglected in this procedure. Us ing these values in Equat ion (3.6), PFsih, t2), and furthermore Re\(t\), w i l l theoretically equal zero, which is an unreasonable solut ion. That p[z(t\), z(t2)] tends to 1 is an unfavorable scenario in evaluating PFs(t\, t2) (Koo and Der Kiureghian 2003; Andr ieu-Renaud et al. 2004). Compl icated by the mult iple dimension of the random variables, the nonlinearity of g, the transformation of the probabil ist ic space, and the dynamic features o f the stochastic process, the evaluation of p[z{t\), z(t2)] has not been clearly and reasonably developed' in the literature, see Der Kiureghian (2000), Ma ie r et al. (2001), and Andr ieu-Renaud et al. (2004) for examples. This paper clearly formulates a concise approximation, which incorporates al l o f the aforementioned factors. 3.3.2 The stochastic linear prediction (SLP) method This method is inspired by the M P F D , or the probabi l ist ical ly geometric centroid o f F. L i and Lence (2007) developed an efficient F O R M - b a s e d approximation to the M P F D in u-space ( u F ) fl(-yg), « F = ^ « (3.17) 48 where uF = the first-order approximation to U F ; <f\( ) - the PDF of the univariate standard normal variable; and a = the negative of the normalized gradient vector of g at u * . Preliminary tests have shown that, in a practical sense, u * can be used in approximating PFS^I, h) for a discrete stationary univariate normal process (Li and Lence 2007). This paper extends the use of u * for approximating PFSOU t2) in generalized discrete cases. FORM divides each u-space at t\ and t2 into a safe domain, S(t\) or S(t2), and a failure domain, F(t\) or F(t2). Figure 3.2 shows a bivariate case. Selecting any point in F(ti), y = (yi, V 2 , . . . , y«)T, the probability of u(fi) falling on the point, P[u(ti) - y], can be written as P[u(ti) = y] = ^(yuy2,...,ya)^u\(ti)Au2(ti)-Aua(ti) (3.18) where #,( ) - the PDF of the n-dimensional standard normal variable; and Au\(ti) = the side length on the t'-th coordinate of a infinitesimal control volume for y. (a) .(b) Figure 3.2. A bivariate example of FORM at two consecutive moments in time (a: at t\ \ b: at t2) With a VARMA model, the value of a stochastic process at a future moment may be predicted based on available data. A common prediction (or forecasting) technique is a linear combination of available data, denoted as SLP in this paper. For details of the technique, see standard time series textbooks (Brockwell and Davis 1996). Assume that the mapping of y in x -space is w . Regardless of its exact form, a one-step SLP based on w , LPi(w), always exists, as long as a time series model for x is well defined. Hence the probability that u(t\) falls on yand LFi ( w ) falls on S(t2) is represented as P[u(t\) = y]-P[LP\(w)eS(t2)]. Cumulatively, once y has ranged over F(t\), one may obtain PFS(H, h) in the following: 49 PFS(^t2)= X P[LPl(wysS(t2M(yl,y2,...yn)Ay ysf(/,) = f P[LPl(w)eS(t2M(yl,y2,...yn)dy l(nP[LPl(w)eS(t2M(yvy2,...yn)dy f . , L>n(>'i.3'2.-3'„)dy =Ef{P[^(w)eS(/2)]}Pf(r,) (3.19) where Ay = AMi(fi)AM2(fi)-"Awn(fi); and Ep{ ) - an expectation over F. It should be noted that P[LPi(w)e S(t2)] is a function of y. Thus Rel(tl)=EF{P[LPl(vr)eS(t2)]} (3.20) Following Li and Lence (2007), Equation (3.20) can be approximated by expanding P[LP\(w)eS(t2)] at u F (fi) , the first-order approximation to the MPFD of u-space at t\. Ignoring second- or higher-order terms, d 'T /?e1(f1)-EF{P[L^(x;(r1))e S(f2)]+ —P[L/»(w)e S(t2)] dy •[y-uF(f,)]} 4 ('i > d T = P{L/>i[(xF(fl)]e S(f2)}+ —P[L/»(w)e S(f2)] dy • E F { [ y - u * F (*,)]} = P{Z^[(xF(f,)]e S(f2)}+ — P[L/>(w)e S(*2)] dy •[EF(y)-uF0*,)] " r ( ' i ) = P{LP1[xl(tl)]eS(t2)} (3.21) where xF (fi) = the mapping of u* (/,) in x-space. It should be noted that E/r(y) = u F (ri). Equation (3.21) is in fact a reliability problem where the random variables are stationary independent standard normal residuals in the time series model of x, denoted as e. The implementation of the SLP method for estimating Rei(t\) is summarized as follows. First, u*(fi) and xF (fi) are obtained; then, g(t2) is rebuilt by replacing x(r2) with LPi[x* (fi)]; and finally, the reliability problem in Equation (3.21) is solved. The physical meaning of Equation (3.21) can be explained in terms of sample paths of u . Originally, PFS^I, h) and Re\(t\) inherently include all sample paths from F{t\) to S(t2); Figure 3.3a illustrates a group of these paths. As shown in Figure 3.3b, Equation (3.21) significantly simplifies the original approach in Figure 3.3a by considering only sample paths from uF(fi) to 50 S(t2). Those sample paths not starting from u F 00 are filtered out. This is a meaningful approximation from a practical viewpoint in that uF(t\) probabilistically represents the entire F(h). Figure 3.3. Sample paths of u inherent in /VsOi, h) and Rei(t\) (a: original; b: filtered) 3.3.3 Summary Both of the BND and SLP methods are applicable to either stationary or non-stationary conditions. Under stationary conditions, the methods are implemented only once to give an estimate of lag-1 resilience that is effective over the entire time horizon. Under non-stationary conditions, however, they must be applied at each time step to delineate the temporal variations of lag-1 resilience, particularly in those cases where seasonal variations are crucial. Both methods require running FORM twice at each time step. The BND method assesses the original performance functions at both of two consecutive moments, while the SLP method evaluates the original performance function at the first moment and a rebuilt one at the later moment. Assume that there are m steps, and that Re\(i) is evaluated for the first m - 1 steps. In this case, the BND method applies FORM m times during the entire procedure, while the SLP 51 method applies 2(m - 1) times. The BND method is thus relatively more computationally efficient than the SLP method. The effectiveness of both innovative methods is case-specific. They may perform differently in different cases, and, thus, it is recommended that one compare both in searching for an appropriate tool for estimating resilience. 3.4 NUMERICAL E X A M P L E 3.4.1 Problem setting: water supply from two sources Two rivers, the White River and the Black River, merge at a community, and provide potential water supply for the community. The water demands of the community vary, and must be satisfied as much as possible by the two rivers. Table 3.1 shows the total monthly water demands of the community, which may include varied components such as industrial, municipal, agricultural, and navigational uses, as well as the monthly mean river flows. It is assumed that each of the White River monthly flows follows a lognormal distribution, while each of the Black River monthly flows follows a normal distribution. The coefficients of variation for these monthly flows are assumed to be 0.1. No long-term trends in water demands and river flows are being considered at this time. Table 3.1. Monthly water demands and mean river flows (Mm3) Month 1* 2 3 4 5 6 7 8 9 10 11 12 Water demand 15 20 20 25 40 45 55 40 25 20 15 12 Mean White 13.6 20 20 21.6 32 38.4 48 37 23 17 14.4 11.2 river flow Black 3.4 5 5 5.4 8 9.6 12 9 6 6 3.6 2.8 (*: Month 1 represents January, and subsequent numbers represent subsequent months.) The monthly flow data of the two rivers may be standardized to random variables with zero mean and unit standard deviation, denoted as s(v, f) = [s\(v, T), s2{y, z)]T, based on: log [ 6 l ( v,.)]-^lo g [ 6 | (: ,r)] } tfUogie, (:,r)]} s 2 ( v , T ) = Q ^ T ) - ^ T ) ] (3.22b) 52 where Q\(v, f) and Q2(v, T) = the original monthly flows of the White and Black Rivers, respectively, at year v and month x, / /( ) = the estimate of the ensemble (or monthly) mean; a( ) = the estimate of the ensemble standard deviation; and s\(v, x) and s2(v, x) = the standardized monthly flows of the White and Black Rivers, respectively, at year v and month x. In this example, it is assumed that the standardized flows are well described by a first-order vector autoregressive model, V A R ( l ) , with constant autoregressive coefficients sm(f + 1): ^m,(r- i - l ) N = A 'sm, (ty + B (eft+iy ysm2(t + 1) ysm2(t)/ ^ 2 ( f + l j , (3.23) where [sm\(t), sm2(t)]T = the vector of standardized total monthly flows of the White and Black rivers generated by the V A R ( l ) model at time t, respectively, or sm(f); A and B = two 2x2 coefficient matrices; and [£\(t + 1), £2(t + 1)]T = a vector of mutually independent standard normal variables at time t + 1, or e(f + 1). Here the time unit is one month. The A and B matrices are determined with the lag-0 and the lag-1 covariance matrices of s, Tsn and Tsi , respectively. For further details regarding the procedure for determining A and B , see standard textbooks (Kottegota 1980; Salas et al. 1980). For illustrative purposes, the covariance matrices defined in Example 4.6 by Kottegota (1980) are assumed here rS() = I 0.845 0.845 1 , Ts, = f0.403 0.397^ 0.272 0.338 In this simple V A R ( l ) , a SLP can be obtained as (Lutkepohl 1991) LP, sm, (t) sm2 (t) = A sm, (t) sm2 (t) + B £, ( ' + 1) K£2(t + l) (3.24) In a given month, a water shortage wil l occur i f the two rivers are unable to meet the water demand in that month. For this water supply problem, reliability, resilience, and vulnerability may be used to probabilistically assess the system performance given monthly changes of system input. In a given month, reliability may be used to indicate the likelihood of a water shortage, vulnerability to indicate the possible magnitude of the shortage, and resilience to indicate the possible length of the shortage. Specifically, lag-1 resilience would indicate the likelihood that, given a water shortage in the previous time period, e.g., the previous month, the two rivers are able to adequately supply the community with water in the current time period, or the current month. (For more details regarding the application of these indices to water supply problems, see Hashimoto et al. (1982)). These probabilistic indicators provide crucial information for 53 determining the capacity of a water supply facility, for planning additional sources such as those from groundwater and reservoir if necessary, and for designing community emergency policies. In practice, this information may be considered along with, for example, water prices for varied water demand components, and water use management options, to evaluate the best alternative for water supply provision. In this example, reliability, Ps, and lag-1 resilience, Re\, of meeting water supply needs of the community are estimated. No long-term trends in water demands and river flows are considered in this example, therefore, Ps and Re\ will not vary annually. Both indicators will vary on a monthly basis, however, due to varying monthly water demands and statistical properties of monthly river flows. Thus the water supply problem is non-stationary. In a given month, t, the water supply problem may be characterized by the following performance function: g(T) = Xl(T) + x2(T)-wd(T) (3.25) where x\(t) and xtit) - the monthly flows of the White and Black Rivers, respectively, at month r, and wd(t) - the water demand at month r. Thus twelve estimates of P a^nd Re\ are obtained for the months of January through December. The monthly flows of the two rivers are considered to be two separate random variables. While it may be more convenient to study the sum of the two river flows as a single random variable, there are sufficient cases in practice that justify separate consideration of both river flows. For example, some water demands within the community might require withdrawals from a specific river, or the water allocation among two source rivers may be necessary and could require consideration of the conditional distributions of both river flows. 3.4.2 Estimates of reliability and resilience The performance function in Equation (3.25) is a linear function of x\ and x2. However, the performance function is a non-linear function of the standardized river flows, because x\ is lognormally distributed and has to be non-linearly transformed. In addition, the joint PDF of x\ and x2 is not available. Hence the PDF of the performance function for this case is difficult to derive analytically. FORM and the BND and SLP methods may thus be considered attractive approaches for estimating Ps and Re\ in this problem and are examined in this example. For comparison sake, MCS is also applied to estimate Ps and Re\. In FORM, design points are obtained with the Hasofer and Lind-Rackwitz and Fiessler (HL-RF) algorithm (Hasofer and Lind 1974; Rackwitz and Fiessler 1978) starting from the u origin, and the convergence criterion is an absolute tolerance of 1CT8. 54 M C S is applied by implementing the following procedure. First, assuming the initial standardized monthly flows to be zero, Equation (3.23) is successively evaluated for 50 time steps to generate 50 warming-up realizations of standardized monthly flows. Following the warming-up realizations, Equation (3.23) continues to generate a 10000-year synthetic data sequence. The data sequence is then transformed back to monthly flows based on Equation (3.22), and these monthly flows are used in Equation (3.25) to generate a synthetic sequence of the performance function. Finally, Ps and Re\ are estimated using the performance function sequence. The variance of the M C S estimates for P$ may be estimated in the following according toMelchers(1999): Var[P s(r)] = [NS(T)—-NS2(T)] (T = 1,...12) (3.26) s NE(NE-\) s NE s where Var( ) = variance; NE - the ensemble size (= 10000 here); and Ns(T) - the number of success realizations at month r. Similarly, the variance of the M C S estimates for Re\ can be roughly estimated as Var[ / te , ( r ) ]« 1 [Nfs(r,T + \)- —l— AT 2 for+ 1)] ( r = l,...12) (3.27) Note that when T - 12, r + 1 wil l return to 1. To compare the efficiencies of the different methods, the number of evaluations of the performance function is used in this analysis as^ -a surrogate for computational cost. Specifically, the computational costs of F O R M and the B N D and SLP methods are primarily incurred in the search for the design points, while those for the M C S are incurred in counting Nf, Ns, and N/s. The B N D method requires the F O R M estimates of P s a n d once these have been estimated, no additional computational costs are incurred. Likewise, once the M C S estimates of Ps have been obtained, no additional computational costs are incurred in M C S for estimating Re\. In contrast, the SLP method not only requires the F O R M estimates of Ps, but additional F O R M evaluations for rebuilt performance functions. 3.4.3 Results Table 3.2 provides all estimates of Ps and Re\ generated for this example. The M C S estimates and their 95% confidence intervals are also shown in Figure 3.4. Table 3.3 lists the number of evaluations of the performance function in each estimation method for each month of the year. 55 Table 3.2. Estimates of Ps and Re\ Month Indicators ; . . 1 2 3 4 5 6 7 8 9 10 11 12 FORM 0891 0.986 0.986 0.772 0.484 0.732 0.802 0.916 0.929 0.917 0.964 0.936 MCS+ 0.894 0.988 0.990 0.784 0.490 0.741 0.809 0.920 0.934 0.923 0.968 0.942 Ps MCS 0.888 0.986 0.987 0.776 0.480 0.732 0.801 0.914 0.929 0.918 0.964 0.938 MCS- 0.882 0.983 0.985 0.768 0.470 0.724 0.793 0.909 0.924 0.913 0.961 0.933 BND 0.949 0.901 0.384 0.272 0.630 0.649 0.810 0.786 0.752 0.874 0.758 0.687 SLP 0.953 0.904 0.382 0.267 0.634 0.653 0.815 0.790 0.756 0.879 0.761 0.690 Re i MCS+ 0.957 0.946 0.523 0.288 0.646 0.667 0.822 0.804 0.763 0.914 0.811 0.718 MCS 0.943 0.896 0.437 0.270 0.633 0.649 0.804 0.776 0.730 0.893 0.768 0.681 MCS- 0.929 0.846 0.350 0.251 0.620 0.631 0.787 0.749 0.698 0.871 0.724 0.644 (MCS+ and MCS-: the upper and lower bounds of 95% confidence intervals for MCS estimates, respectively.) 0.2 I ' ' ' ' < 1 1 ' 1 1 ; 1 2 3 4 5 6 7 8 9 10 11 12 Month Figure 3.4. MCS estimates and their 95% confidence intervals (The bars indicate the confidence intervals.) Table 3.3. Number of performance function evaluations Methods 1 2 3 4 Month 5 6 7 8 9 10 11 12 FORM 5 6 6 5 4 5 5 5 5 5 5 5 MCS 10000 for each month SLP* 5 5 5 5 5 5 5 5 5 5 5" 5 (*: The additional cost for estimating Re\.) 56 Given the range of the 95% confidence limits for the MCS estimates of Ps, the ensemble size of 10000 is considered to be sufficiently large to control sampling errors of these MCS estimates. As shown in Table 3.2, the FORM estimates of Ps are always within the 95% confidence limit of the MCS estimates of Ps, indicating that FORM is rather effective for this example. The BND and SLP methods produce very close estimates of Re\, and agree with those generated by MCS, falling within the 95% confidence interval in all cases. Deeper investigations may be needed to reveal their relationships. The scatter of the MCS estimates of Re\, which is generally greater than that for the MCS estimates of Ps, reflects the efficiency of the sampling method in estimating resilience. As shown in Figure 3.4, the range of the 95% confidence limits generally increases as Ps increases. For example, in May Ps is relatively low, and the range of the 95% confidence limits of Re\ is 0.026, or the coefficient of variation (COV) is 4%. In March, on the other hand, with a relatively high value of Ps, this range is 0.173, or the COV is 40%. The phenomenon may be explained by examining Equation (3.3). When Ps(r) is high, the corresponding value of NJ(T) will be small. Hence the accuracy of Equation (3.3) decreases, and more generations are needed to enhance accuracy, reducing the efficiency of MCS. The results indicate that a high reliability does not necessarily imply a high resilience. For example, the water supply reliability in March is quite high due to medium water demands, and adequate water availability at this time. If a water shortage occurs in March, however, it will be difficult to meet the subsequent water demands of April, and the shortfall will likely last two months. The lag-1 resilience depends on the system performance at two consecutive periods, rather than simply the likelihood of system functioning at any single moment. It may be used by system operators to understand system performance, and to make informed decisions regarding the allocation of shortfalls over time considering the trade-off between meeting demands today and meeting them in the near future. As shown in Table 3.3, the BND and SLP methods are more computationally efficient than MCS. The HL-RF algorithm converges quickly in searching for design points for each month. In estimating Re\ throughout the year, the BND and SLP methods evaluate the performance function 61 times and 121 times, respectively, while MCS evaluates performance function 120000 times. Considering their accuracy compared with MCS, in this example, the BND and SLP methods may be used as an efficient tool for estimating Re\. The system performance may be improved by adjusting water demand or developing additional supply, depending on the preferences of the community. A change in demand or supply at any time t influences Ps at t, and Re\ at both t and t - 1. As an example of the 57 usefulness of estimates of Rei and Ps for selecting among various options, consider the case in which the system operator wishes to improve the system performance in Apr i l . It is assumed that an additional water supply of 2 M m may be provided to the community from another source in either April or May. Two alternatives are to be considered. In Alternative 1, the additional water supply is utilized in Apri l , and in Alternative 2 it is utilized in May. F O R M and the B N D method are applied to estimate Ps and Re\ under these alternatives, and Table 3.4 provides these results. Since the additional water supply does not influence system performance in other months, only the results for March, Apri l , and May are provided. The estimates indicate that utilizing the additional source in a given month wil l improve Ps in that month and Rei in the previous month, but decrease Rei in that month. In Alternative 1, for instance, with an additional supply in Apr i l , the likelihood of meeting the water demand in Apr i l increases from 0.772 to 0.943. Meanwhile, given a water shortage in March, the likelihood of no shortage in Apri l increases from 0.384 to 0.730. However, given a water shortage in Apri l , the likelihood of no shortage in May decreases from 0.272 to 0.182. This is an understandable result. With the additional supply and the resultant high reliability, a water shortage in Apri l might result simply from an extreme low-flow event in Apr i l . Following the extreme event, the two rivers may hardly provide adequate supply in May, due to the lag-1 covariance structure of the river flows. Table 3.4. Impacts of an additional water supply Indicators Water supply schemes 3 Month 4 5 Original 0.986 0.772 0.484 Alternative 1 0.986 0.943 0.484 Alternative 2 0.986 0.772 0.686 Ret Original Alternative 1 Alternative 2 0.384 0.730 0.384 0.272 0.182 0.482 0.630 0.630 0.572 3.5 C O N C L U S I O N S Resilience reflects the recovery capacity of a system that is aided by natural hydrological or meteorological cycles. The recovery capacity is, thus, important in those systems largely influenced by cyclic natural conditions, such as aquatic environmental or ecological systems. For more information regarding the importance of resilience, see Maier et al. (2002). 58 This paper has proposed two practical approximation methods for lag-1 resilience, based on FORM and the periodic VARMA model. They are designed for accommodating complex scenarios, i.e., multivariate, non-stationary, non-linear, and non-normal problems. Although developed in the discrete environment, these two methods may be also applicable to continuous cases, because continuous stochastic processes can be generally discretized to discrete stochastic processes. See Andrieu-Renaud et al. (2004) for details regarding the discretization procedure. In this work the effectiveness of these methods has been validated for a multivariate, non-stationary, non-linear, and non-normal problem. However, this effectiveness is conditional; that is, the accuracy of these methods will be reduced if FORM is not effective for a given problem, and their feasibility may be affected if a VARMA model is not available due to inadequate data. Both the BND and the SLP method make substantial contributions to the solution of problems regarding resilience. The former theoretically formulates an approximation, while the latter creates an innovative application of the MPFD. Both methods are substantially more computationally efficient than sampling methods, and may provide competitive options for estimating resilience in cases where FORM is effective and a periodic VARMA model is available. Moreover, these two practical methods may advance the inclusion of resilience as a goal in risk-based water resources management programs, especially for situations where seasonal effects are a major concern. 59 R E F E R E N C E S Andrieu-Renaud, C , Sudret, B., and Lemaire, M. (2004). "The PHI2 method: a way to compute time-variant reliability." Reliability Engineering & System Safety, 84(1), 75-86. Ang, A. H. S., and Tang, W. H. (1984). Probability Concepts in Engineering Planning and Design, Vol. II: Decision, Risk, and Reliability, John Wiley, New York. ASCE Task Committee on Sustainability Criteria. (1998). Sustainability Criteria for Water Resource Systems, American Society of Civil Engineers. Bras, R. L., and Rodriguez-Iturbe, I. (1985). Random Functions and Hydrology, Addison-Wesley Publishing Company, Reading, MA. Brockwell, P. J., and Davis, R. A. (1996). Introduction to Time Series and Forecasting, Springer-Verlag, New York, NY. Chiang, P.-C, and Gates, T. K. (2004). "Strategic River Water Quality Planning Using Calibrated Stochastic Simulation." ASCE Journal of Water Resources Planning and Management, 130(3), 215-231. Coles, S., and Tawn, J. (2005). "Seasonal Effects of Extreme Surges." Stochastic Environmental Research and Risk Assessment, 19(6), 417-427. Deo, M. C , Sherief, M., and Sarkar, A. (1997). "Wave Height Estimation Using Disaggregation Models." ASCE Journal of Waterway, Port, Coastal, and Ocean Engineering, 123(2), 63-67. Der Kiureghian, A. (2000). "The Geometry of Random Vibrations and Solutions by FORM and SORM." Probabilistic Engineering Mechanics, 15(1), 81-90. Durrans, S. R., Eiffe, M. A., Thomas Jr., W. O., and Goranflo, H. M. (2003). "Joint Seasonal/Annual Flood Frequency Analysis." ASCE Journal of Hydrological Engineering, 8(4), 181-189. Fiering, M. B. (1976). "The Role of Systems Analysis in Water Program Development." Natural Resources Journal, 16, 759-771. Fiering, M. B. (1982). "Alternative Indices of Resilience." Water Resources Research, 18(1), 33-39. Fowler, H. J., Kilsby, C. G., and O'Connell, P. E. (2003). "Modeling the Impacts of Climate Change and Variability on The Reliability, Resilience, and Vulnerability of A Water Resource System." Water Resources Research, 39(8), 1222. Haan, C T . (2002). Statistical Methods in Hydrology, Iowa State Press, Ames. Harr, M. E. (1987). Reliability-Based Design in Civil Engineering, McGraw-Hill Book Company, New York. 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. C. (1974). "Exact and Invariant Second-Moment Code Format." ASCE Journal of the Engineering Mechanics Division, 100(1), 111-121. Holling, C. S. (1973). "Resilience and Stability of Ecological Systems." Annual Review of Ecological Systems, 4, 1-23. Kac, M. (1943). "On The Average Number of Real Roots of A Random Algebraic Equation." Bulletin of the American Mathematical Society, 49, 314-320. Karamouz, M., Szidarovszky, F., and Zahraie, B. (2003). Water Resources Systems Analysis, CRC Press LLC, Boca Raton, FL. Kjeldsen, T. R., and Rosbjerg, D. (2004). "Choice of Reliability, Resilience and Vulnerability Estimators for Risk Assessments of Water Resources Systems." Hydrological Sciences Journal, 49(5), 755-767. 60 Koo, H., and Der Kiureghian, A. (2003). "FORM, SORM, and Simulation Techniques for Nonlinear Random Vibrations." UCB/SEMM-2003/01, Department of Civil and Environmental Engineering, University of California, Berkeley, California. Kottegota, N. T. (1980). Stochastic Water Resources Technology, The Macmillan Press Ltd., Hong Kong. Li, C. Q., Lawanwisut, W., and Zheng, J. J. (2005). "Time-Dependent Reliability Method to Assess the Serviceability of Corrosion-Affected Concrete Structures." ASCE Journal of Structural Engineering, 131(11), 1674-1680. Li, Y., and Lence, B. J. (2006). "On the Applicability of Rice's Formula in Stochastic Hydrological Modeling." ASCE Journal of Hydrological Engineering, in review. Li, Y., and Lence, B. J. (2007). "Mean Points of Failure and Safety Domains: Estimation and Application." Probabilistic Engineering Mechanics, 22(3), 257-266. Loucks, D. P. (1997). "Quantifying Trends in System Sustainability." Hydrological Sciences Journal, 42(4), 513-530. Liitkepohl, H. (1991). Introduction to Multiple Time Series Analysis, Springer-Verlag, Berlin. Luz, L. D., and Loucks, D. P. (2003). "Developing Habitat Suitability Criteria for Water Management: A Case Study." International Journal of River Basin Management, 1(4), 283-295. Madsen, H. O., Krenk, S., and Lind, N. C. (1986). Methods of Structural Safety, Prentice-Hall, Englewood Cliffs, NJ. Maier, H. R., Lence, B. J., and Tolson, B. A. (2002). "The Role of Reliability, Vulnerability and Resilience in the Management of Water Quality." the 27th Hydrology and Water Resources Symposium, The Institution of Engineers, Australia, Melbourne, Australia. Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (2001). "First-Order Reliability Method for Estimating Reliability, Vulnerability and Resiliency." Water Resources Research, 37(3), 779-790. Melchers, R. E. (1999). Structural Reliability Analysis and Prediction, John Wiley & Sons Ltd., Chichester. Nordin, C. F., and Rosbjerg, D. (1970). "Applications of Crossing Theory in Hydrology." Bulletin of the International Association of Scientific Hydrology, 15(1), 27-43. Putcha, C. S., and Patev, R. C. (1998). "Time-Variant Reliability Assessment." Uncertainty Modeling and Analysis in Civil Engineering, B. M. Ayyub, ed., CRC Press LLC, Boca Raton, FL, 73-97. Rackwitz, R., and Fiessler, B. (1978). "Structural Reliability under Combined Random Load Sequences." Computers & Structures, 9, 489-494. Rice, S. O. (1944). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 23(3), 282-332. Rice, S. O. (1945). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 24(1), 46-156. Rodnguez-Iturbe, I. (1968). "A Modern Statistical Study of Monthly Levels of The Orinoco River." Bulletin of the International Association of Scientific Hydrology, 13(4), 25-40. Rodriguez-Iturbe, I. (1969). "Applications of The Theory of Runs to Hydrology." Water Resources Research, 5(6), 1422-1426. Rosbjerg, D. (1977). "Crossing and Extremes in Dependent Annual Series." Nordic Hydrology, 8(5), 257-266. Salas, J. D. (1996). "Analysis and Modeling of Hydrologic Time Series." Water Resources Handbook, L. W. Mays, ed., McGraw-Hill, 19.01-19.72. 61 Salas, J. D., Chung, C.-H., and Cancelliere, A. (2005). "Correlations and Crossing Rates of Periodic-Stochastic Hydrologic Processes." ASCE Journal of Hydrological Engineering, 10(4), 278-287. Salas, J. D., Chung, C.-H., and Fernandez, B. (2001). "Relating Autocorrelations and Crossing Rates of Continuous- and Discrete-Valued Hydrologic Processes." ASCE Journal of Hydrological Engineering, 6(2), 109-118. Salas, J. D., Delleur, J. W., Yevjevich, V., and Lane, W. L. (1980). Applied Modeling of Hydrological Time Series, Water Resources Publications, Littleton, CO. Sen, Z. (1976). "Wet and Dry Periods of Annual Flow Series." ASCE Journal of the Hydraulics Division, 102(HY10), 1503-1514. Shooman, M. L. (1968). Probabilistic Reliability: An Engineering Approach, McGraw-Hill Book Co., New York. Sundararajan, C. (1995). "Probabilistic Structural Mechanics Handbook." Chapman & Hall, New York. Tick, L. J., and Shaman, P. (1966). "Sampling Rates and Appearance of Stationary Gaussian Processes." Technometrics, 8(1), 91-106. Tung, Y. K., Yen, B. C , and Melching, C. S. (2006). Hydrosystems Engineering Reliability Assessment and Risk Analysis, McGraw-Hill, New York. UNESCO/IHP. (1999). Sustainability Criteria for Water Resource Systems, Cambridge University Press, Cambridge, UK. Weiler, K, Walter, M. T., Walter, M. F., Brooks, E. S., and Scott, C. A. (2000). "Seasonal Risk Analysis for Floodplains in the Delaware River Basin." ASCE Journal of Water Resources Planning and Management, 126(5), 320-329. 62 C H A P T E R 4 O N T H E APPLICABIL ITY O F RICE'S F O R M U L A IN S T O C H A S T I C H Y D R O L O G I C A L M O D E L I N G A version of this chapter has been submitted to ASCE Journal of Hydrological Engineering. The paper is titled as On the Applicability of Rice's Formula in Stochastic Hydrological Modeling by Y. Li and B. J. Lence. 63 P R E F A C E Rice's formula is used to measure the instantaneous crossing rate for a given threshold level. The formula was originally derived for a continuous process. The formula essentially utilizes a stochastic linear prediction (SLP) of a stochastic process. Rice's formula has been used in varied engineering problems, such as telecommunications, structural dynamics, signal processes, and material fatigues. The evaluation of the crossing probability for a given threshold level of a discrete process is an important modeling issue in stochastic hydrology. The application of Rice's formula for this issue is also available in the hydrological literature. However, hydrological modeling is typified by various discretized representations of stochastic processes. Whether the formula is truly suitable for estimating two-state lag-1 transition probabilities and the associated resilience in hydrological modeling is the theme of this chapter. This chapter examines in depth the discrete application of Rice's formula to estimating the up-crossing probability in typical discrete stationary univariate Gaussian processes. The results based on Rice's formula are shown to follow a biased pattern relative to exact solutions, and to perform poorly in several important scenarios, particularly for estimating resilience in these numerical experiments. The results of this chapter indicate that Rice's formula may be an inappropriate methodology for approximating the two-state lag-1 transition probabilities in discrete processes, as well as the associated resilience. The reasonableness of the two approximation methods for resilience presented in Chapter 2 is thus indirectly justified. 64 4.1 INTRODUCTION The classical Rice's formula was originally developed to estimate the instantaneous crossing rate of a continuous stationary univariate Gaussian process beyond a given threshold level (Kae 1943; Rice 1944, 1945). As a significant contribution to the theory of crossings or runs, the formula has been widely applied for estimating crossing probability for a given threshold and other related indicators in various engineering problems, such as telecommunications, structural dynamics, signal processes, and material fatigues. In stochastic hydrology, it has been employed in modeling the time series of a single hydrological input, e.g. a flow rate or water level (Rodriguez-Iturbe 1968, 1969; Nordin and Rosbjerg 1970; Bras and Rodriguez-Iturbe 1985; Salas et al. 2001; Salas et al. 2005). Hydrologists have also warned that Rice's formula, as a limit expression for continuous cases, may be an inappropriate tool to describe hydrological processes that are primarily comprised of discrete observations in practice (Nordin and Rosbjerg 1970; Kottegota 1980). This technical note demonstrates the shortcomings of this formula in estimating crossing probability and resilience for typical discrete processes. 4.2 RICE'S FORMULA Expressed in various formats, Rice's formula can be written as follows: - ( C - f t t ) 2 n -(C-Mx f y + =—^JLe W = 1 r ~ R X " ( Q ) - | O . 5 c 2 [ ^ ( 0 ) - ^ 2 ] ,4 j s ° 27tOx 271 Rx(0)-{ix2 where f = time; X = X(t), a continuous stationary normal process; C = a fixed threshold level; vc = the instantaneous up-crossing rate of X over C, a constant due to the stationarity of X; X = X (t), the derivative process of X; jUx = the mean of X; Rx = the autocorrelation function of X; Rx" = the second-order derivative function of Rx', <Jx = the standard deviation of X; and ox = the standard deviation of X . The derivation of the formula is not hard to understand, see Rice (1944; 1945) or Melchers (1999) for details. Note that down-crossings and up-crossings are studied in the same way, and that this technical note focuses on up-crossings without any loss of generality. Obviously, the mean number of up-crossings during a given period should be obtained by integrating vc+ over the period. If the period is small, the following expression will hold at any moment fo: Nu(t0, f0 + St) = Pu(t0, f0 + St) ~ vc+St (4.2) 65 where St = an infinitesimal time interval; Nu(to, to + dt) - the mean number of up-crossings during [to, to + St]; and Pu(to, to + St) = the probability of an up-crossing occurring during [to, to + St]. Also see Koo and Der Kiureghian (2003) for the derivation of Equation (4.2). Two important conditions must be noted regarding Rice's formula. First, it is a limit in the continuous time domain, thus Equation (4.2) is applicable only when St is adequately small. Second, Equation (4.1) requires differentiating Rx, but this task may encounter difficulties in reality. As an example, Figure 4.1 illustrates three scenarios of Rx. In scenario (a), Rx"(0) can be obtained analytically, or approximated numerically in a reasonable manner. In scenario (b) or (c), however, Rx is non-differentiable at t - 0, thus no analytical solution of /?A-"(0) exists, and a numerical approximation is not realistic. If the two conditions are not met, Equation (4.2) may cause serious errors in modeling discrete processes, as will be illustrated in the next section. In addition, X , and Rx ' are hard to obtain, and even hard to understand for average engineers. Rice's formula and associated methods may thus be difficult to implement in practice, as also pointed out by Putcha and Patev (1998). See Li et al. (2005) for the complicated implementation of Rice's formula. (a) (b) (c) Figure 4.1. Examples of Rx(t) 4.3 C R O S S I N G S I N D I S C R E T E P R O C E S S E S Most hydrological data are recorded or used in a discrete format, though they are intended to describe continuous phenomena. Rain gauges record rainfall intensity usually in five-minute intervals, and stream flow rates are widely represented by average values (annually, monthly, biweekly, weekly, or daily). In such cases, lag-1 crossing probabilities for stationary univariate normal processes can be estimated with a joint bivariate distribution (Tick and Shaman 1966) PUl(t0) = £ £°&[x(t0),x(t0 + l)^(l)]dx(/0)dx(r0 + 1) (4.3) 66 where x(to) = the value of X at to; x(to + 1) = the value of X at one time interval after to; p(l) = the lag-1 autocorrelation coefficient of X; Pu\{to) = the probability of an up-crossing occurring during the time period [to, to + 1], or the lag-1 up-crossing probability at to; and ) - the probability density function (PDF) of a bivariate normal distribution. Because of the discreteness of the X process, the stochastic model of X is also built on a discrete basis. As a result, Equation (4.3) is the exact solution for PwiOo) in the discrete time domain, where the time unit is case-specific, depending on modeling requirements. Despite the existence of Equation (4.3), Rice's formula is still applied in stochastic hydrological modeling (Rodriguez-Iturbe 1968, 1969; Bras and Rodriguez-Iturbe 1985; Salas et al. 2001; Salas et al. 2005). In a logical sense, however, such applications of Rice's formula may be unnecessary and unreasonable. In these applications, Equation (4.2) is generally simplified to Pu[(t0) = vc+-(t0+l-tQ) = vc+ (4.4) A key to applying Rice's formula relies on the estimation of ax and Px"(0). Assuming Rx to be twice differentiable at t = 0, /?x"(0) can be numerically approximated with a central difference scheme (Rodriguez-Iturbe 1968) RX"(0) = 2RX(0)-2RX(\) (4.5) . However, the reasonableness of this approximation may not been proven in some cases. Likewise, a backward difference scheme can be employed for X (t) as follows: X(t0)= * ' o ) - * < b - l ) = X ( , Q } _ x ( ? o _ 1 } ( 4 - 6 ) Oo 1) Consequently, cr.2 = <r2[x(t0)] + a2[x(tQ -1)] + 2 cov[jc(r0), -x(t0 -1)] = 2o2-2a2pt\) (4.7) where cov( ) = autocovariance. Both approximations of /?x"(0) and ax lead to the following same result: Pul(tQ) = ^ 0-e£^' (4.8) 4.4 NUMERICAL EXAMPLES Numerical examples are herein provided to compare the applications of Equations (4.8) and (4.3). In addition, estimates of lag-1 resilience based on Equations (4.8) and (4.3) are provided 67 for a deeper comparison. The lag-1 resilience at t - to, Re\(t0), represents the following conditional probability: Re](t0) = P[x(t0 + l)>C/x(t0)<C] _P[x(tQ)<C,x(t0 + l)>C] PM,(f 0) P[x(t0)<C] (4.9) P[x(t0) < C] where P[ ] = probability. The resilience measures the recovery capacity of repairable risk systems, such as aquatic environmental systems or water supply systems, from unsatisfactory states to satisfactory states (Hashimoto et al. 1982; Maier et al. 2001; Fowler et al. 2003). Figure 4.2 shows a stationary normal process X with zero mean and unit standard deviation. Equations (4.8) and (4.3) and the resulting resilience are tested for different values of C and p{\). In this stationary example, either the lag-1 up-crossing probability or the lag-1 resilience is constant over time. A the PD ^ normal F of standard distribution | y \ to I ro+1 \ \ Figure 4.2. A stationary Gaussian process X Figure 4.3 illustrates the results of Pui(t0). Some parts of the results are also available in Nordin and Rosbjerg (1970), Salas et al. (2001), and Li and Lence (2007). Rice's formula converges to the exact solution as p{\) approaches 0.9. However, Rice's formula follows a pattern that is slightly biased from the exact solution based on Equation (4.3). The difference between Equations (4.3) and (4.8) generally increases as p(\) decreases, particularly when C equals zero. In addition, Rice's formula provides a poor estimate when both C and /XI) are zero, 68 a very representative scenario that must be tested in building stochastic models. This important discrepancy can also be observed in Nordin and Rosbjerg (1970) and Li and Lence (2007). 0.45 s a . 0.15 Figure 4.3. Computational results for Pu\(to) Figure 4.4 illustrates the results of Re\(to). The solution based on Rice's formula again converges to the analytical solution as p(\) approaches 0.9. The errors caused by Rice's formula increase as p(l) decreases. The errors are typically greater when C becomes negative. Rice's formula is unable to accurately estimate Rei(t0) when both C and p{\) are zero. A serious problem in applying Rice's formula is observed for cases where C < 0 in that some of its results are larger than 1.0. Such results are not reasonable, for resilience represents a probabilistic content which must not exceed 1.0. -0.3 0 M D Figure 4.4. Computational results for Re [(to) 69 4.5 CONCLUSIONS The discrete solutions of Rice's formula and the analytical solutions to the lag-1 up-crossing probability and the lag-1 resilience in discrete processes have been thoroughly examined using simple but typical numerical examples. It is found that Rice's formula cannot give satisfactory results in an overall sense. Particularly, serious errors may occur in estimating resilience. The findings in this note warn hydrologists that, in general, Rice's formula may not be a proper tool in discrete stochastic hydrological modeling. The warning, however, is not to omit successful applications of the formula in some cases. The two conditions mentioned earlier might be the major cause for the inapplicability. The non-differentiability of Rx may reduce or eliminate the applicability of Rice's formula even in continuous processes. Equation (4.3) may be employed without hesitance for estimating the up-crossing probability in a discrete stationary univariate process, though a fine time unit could improve the performance of Rice's formula. To the authors' knowledge, Rice's formula has not been used for estimating resilience in water resources systems. The findings in this note provide an argument against employing Rice's formula for estimating resilience in discrete processes. Furthermore, the findings help in selecting the proper method for estimating transition probabilities and resilience in non-stationary, non-linear and multivariate hydrological processes, or other potential processes in engineering practice. 70 R E F E R E N C E S Bras, R. L., and Rodriguez-Iturbe, I. (1985). Random Functions and Hydrology, Addison-Wesley Publishing Company, Reading, MA. Fowler, H. J., Kilsby, C. G., and O'Connell, P. E. (2003). "Modeling the Impacts of Climate Change and Variability on The Reliability, Resilience, and Vulnerability of A Water Resource System." Water Resources Research, 39(8), 1222. 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. Kac, M. (1943). "On The Average Number of Real Roots of A Random Algebraic Equation." Bulletin of the American Mathematical Society, 49, 314-320. Koo, H., and Der Kiureghian, A. (2003). "FORM, SORM, and Simulation Techniques for Nonlinear Random Vibrations." UCB/SEMM-2003/01, Department of Civil and Environmental Engineering, University of California, Berkeley, California. Kottegota, N. T. (1980). Stochastic Water Resources Technology, The Macmillan Press Ltd., Hong Kong. Li, C. Q., Lawanwisut, W., and Zheng, J. J. (2005). "Time-Dependent Reliability Method to Assess the Serviceability of Corrosion-Affected Concrete Structures." ASCE Journal of Structural Engineering, 131(11), 1674-1680. Li, Y., and Lence, B. J. (2007). "Mean Points of Failure and Safety Domains: Estimation and Application." Probabilistic Engineering Mechanics, 22(3), 257-266. Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (2001). "First-Order Reliability Method for Estimating Reliability, Vulnerability and Resiliency." Water Resources Research, 37(3), 779-790. Melchers, R. E. (1999). Structural Reliability Analysis and Prediction, John Wiley & Sons Ltd., Chichester. Nordin, C. F., and Rosbjerg, D. (1970). "Applications of Crossing Theory in Hydrology." Bulletin of the International Association of Scientific Hydrology, 15(1), 27-43. Putcha, C. S., and Patev, R. C. (1998). "Time-Variant Reliability Assessment." Uncertainty Modeling and Analysis in Civil Engineering, B. M. Ayyub, ed., CRC Press LLC, Boca Raton, FL, 73-97. Rice, S. O. (1944). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 23(3), 282-332. Rice, S. O. (1945). "Mathematical Analysis of Random Noise." Bell System Technical Journal, 24(1), 46-156. Rodriguez-Iturbe, I. (1968). "A Modern Statistical Study of Monthly Levels of The Orinoco River." Bulletin of the International Association of Scientific Hydrology, 13(4), 25-40. Rodriguez-Iturbe, I. (1969). "Applications of The Theory of Runs to Hydrology." Water Resources Research, 5(6), 1422-1426. Salas, J. D., Chung, C.-H., and Cancelliere, A. (2005). "Correlations and Crossing Rates of Periodic-Stochastic Hydrologic Processes." ASCE Journal of Hydrological Engineering, 10(4), 278-287. Salas, J. D., Chung, C.-H., and Fernandez, B. (2001). "Relating Autocorrelations and Crossing Rates of Continuous- and Discrete-Valued Hydrologic Processes." ASCE Journal of Hydrological Engineering, 6(2), 109-118. Tick, L. J., and Shaman, P. (1966). "Sampling Rates and Appearance of Stationary Gaussian Processes." Technometrics, 8(1), 91-106. 71 C H A P T E R 5 CONCLUSIONS This thesis is an engineering probabilistic study that is applicable to risk analysis for hydrotechnical systems. The prime focus is the estimation and analysis of two-state lag-1 transition probabilities and associated resilience in discrete processes, which are common in an array of hydrotechnical problems, such as water allocation and quality management. The methods developed here are also applicable to other civil engineering fields, such as structural reliability analysis under dynamic conditions. 5.1 SUMMARY The major contributions of this dissertation are summarized here. A systematic summary of the existing techniques for estimating two-state lag-1 transition probabilities in discrete processes The evaluation of the two-state lag-1 transition probabilities and derivatives is a long-time research topic in the analysis and application of telecommunications, statistics, signal process, structural reliability, hydrological modeling, and many other disciplines. For decades, numerous works have been published for estimating Pis related to the transition probabilities. These efforts are interrelated, because the Pis stem from the same root. On the other hand, different estimation methodologies have been applied due to the varied features of different problems. The features may be reflected in various ways, and may be identified, for example, by addressing the following questions. Is the stochastic process to be studied discrete or continuous? Does the process involve single or multiple random variables (or component processes)? Is the process a linear or non-linear function of component processes? Are component processes normal or non-normal? Are component processes inter-correlated or not? Are component processes auto-correlated or not? In general, nevertheless, no profound and broad overview regarding these methodologies across different engineering disciplines has been provided in the literature. Meanwhile, the degree of technical advances in this topic for the past several decades has been varied over different disciplines, being dramatic in some disciplines, and nominal in others. As a result, from a modern viewpoint, these research efforts appear somewhat isolated and narrowly focused, in terms of the methodologies for estimating the transition probabilities. When faced with so many interrelated and varied methods, one may get lost, and may ask why there are so many ways, how different they are in nature, and which way we should go for a specific problem. 73 To improve this tangled situation, and to clarify the relationship among different methodologies, this dissertation provides a broad review on the methods for estimating the two-state lag-1 transition probabilities in discrete processes. The review basically reflects the historical development and the contemporary status of this topic, and comments on practical methods for estimating transition probabilities and resultant resilience considering practicality and generality. The review categorizes the existing techniques, and especially identifies two major approaches, which have been widely applied in the literature, i.e. the SLP-based method in conjunction with Rice's formula, and the BND-based method. Their shortcomings are also addressed. The review not only constructs the foundation for the original theoretical development in this research, but also, more importantly, provides an informative setting for future research in similar topics. The essential points of the review are separately presented in this chapter, and substantial details have been presented in the previous chapters. Identification of a new fundamental property of reliability problems, Le. MPFD and MPSD The load-resistance interference mechanism for reliability problems has been intensively studied for decades. With the introduction of modern reliability techniques such as FORM and SORM over the past three decades, the geometric analyses of the interference mechanism have brought substantial advances in engineering probabilistic study. This thesis makes a significant contribution to reliability analysis by introducing MPFD and MPSD. (See Chapter 2 for details.) Both points represent the geometric centroid of the standardized probabilistic space (u). A practical approximation for them, Equations (2.7) and (2.9), is developed herein based on FORM, making the application of both points feasible in practice. These points provide a new way of thinking about reliability problems, and may open new research fields, such as an innovative SLP-based method for estimating the two-state transition probabilities that is developed herein. More topics on the application of the mean points are expected to be identified, and are discussed in Section 5.3. It is worth stressing again that the development of the approximation and the innovative SLP-based method largely benefit from a modern reliability technique, FORM, which is most intensively employed in structural reliability analysis. Development of a new SLP-based method in conjunction with MPFD and MPSD As mentioned above, an innovative SLP-based method for estimating the two-state transition probabilities is developed herein on the basis of MPFD and MPSD. The method is the first practical application of the mean points. (See Chapters 2 and 3 for details.) 74 This is a reasonable application of the mean points, and the physical implication of the method is easy to understand. The stochastic process of performance function originally consists of numerous sample paths. The nature of two-state transition probabilities is the portion of sample paths of the stochastic process that travel from one system state at an earlier time moment into the other system state at a later time moment. The sample paths from one domain may be approximately represented by those starting from the mean point of that domain, in a probabilistic sense. By this philosophy, the new SLP-based method investigates only those sample paths, significantly simplifying the analysis of the transition probabilities. This achievement reflects a combination of techniques in different disciplines. The basis of the method comprises two parts, i.e. FORM and time series modeling. The former was developed mainly by practitioners in structural engineering, while the latter is a common tool in stochastic hydrology. Development of a practical approximation approach for estimating the lag-1 autocorrelation coefficient of g process Again based on FORM and time series modeling, this thesis theoretically formulates a concise approximation approach for the autocorrelation coefficient. (See Equation (3.14), Chapter 3 for details.) The autocorrelation coefficient is a key in the BND-based method for estimating two-state lag-1 transition probabilities. The approximation approach systematically and explicitly incorporates the lag-1 covariance structure of random variables. It is straightforward to implement, unlike existing approaches in which the autocorrelation coefficient may be only descriptively expressed, for example, those using the KL representation method for estimating the autocorrelation coefficient. This research thus improves the current BND techniques. Investigation of the application of Rice's formula As mentioned before, there are different methodologies for estimating the two-state transition probabilities. One of them is the SLP-based method used in combination with Rice's formula, which was introduced in the 1940's. It is theoretically different from another common methodology, the BND-based method. The BND-based method has a clearer and more understandable principle than the Rice's formula-based method, and has been widely used since the 1970's. However, both methods are employed in parallel by structural engineering researchers. In addition, the author found that Rice's formula may produce unsatisfactory 75 modeling results for discrete processes, even in very simple cases. These results in fact have been reported for almost 40 years in the literature, but no deep studies have been conducted to examine the discrete application of the formula, though the formula is still used in discrete processes today. To investigate the applicability of Rice's formula in discrete processes, this thesis conducts a series of numerical experiments based on a typical stochastic modeling problem. (See Chapter 4 for details.) The experiments fully demonstrate the modeling applicability and limitations of the formula, and indicate that Rice's formula may not be robust in discrete cases. The investigation may have an important influence on the discrete application of Rice's.formula: Development of a risk analysis approach for a given type of water supply problem To demonstrate the proposed resilience methods, Chapter 3 (Section 3.4) presents a hypothetical case study of the management of a local water supply. The example may be instructive for cities, municipalities, or regions with water supply sources from multiple streams. Most rivers have an obvious annual cycle, and evaluating typical Pis on a seasonal basis will help decision-makers to understand the water shortage risks of such a system, and to identify robust mitigation strategies, as demonstrated in the example. In the summer of 2006, Chongqing, China experienced an unprecedented drought. Here, the Yangtze River and one of its tributaries, the Jialing River, are major water supply sources to municipal, agricultural, and industrial demands. The demands are normally satisfied with water withdrawn from one of the two rivers. Usually, summer is the high-flow season for both rivers. In the summer of 2006, however, both rivers simultaneously experienced extreme low flows, which lasted for about two months, and the ensuing drought resulted in a serious economical loss to the city. The example presented in Chapter 3 may provide an appropriate approach for analyzing the water supply safety of the City of Chongqing, and of other municipalities subjected to similar water supply conditions. 5.2 CONCLUSIONS The following conclusions can be drawn based on the contributions of this thesis. • The two-state lag-1 transition probabilities in discrete nonlinear multivariate processes may be effectively estimated within a framework of a time series model. Point process models are not suitable in this case. 76 • Currently, there are two major methodologies for estimating the transition probabilities, 1) the BND-based method, and 2) the SLP-based method in conjunction with Rice's formula. • In the existing BND-based method, the techniques employed to estimate the lag-1 autocorrelation coefficient of g process may significantly increase the number of random variables. • The current Rice's formula-based method may cause considerable errors in modeling discrete processes, and its implementation procedure is complicated. • The MPFD and MPSD are among the fundamental properties of reliability problems, and may be useful in practical applications. • The MPFD and MPSD may be easily approximated using FORM. • The MPFD and MPSD can be used to build a SLP for estimating two-state lag-1 transition probabilities in discrete processes, as well as resilience. • The lag-1 autocorrelation coefficient of the g process can be conveniently approximated based on FORM, without increasing random variables. The approximation leads to a new BND-based method for estimating the transition probabilities and the resultant resilience. • The seasonal evaluation of Pis can provide realistic descriptions of risks for a municipal water supply system with two or more rivers as the main water source. 5.3 FUTURE WORK More applications of MPFD and MPSD will be examined in the future. For example, as mentioned in Chapter 2, the MPFD can be a good candidate for the center of sampling density function in an importance sampling scheme. Usually, the design point is used as the mean of the sampling density function in importance sampling. As the MPFD is located inside the failure domain, a sampling density function centered on the MPFD might sample more realizations in the failure domain than one centered on the design point. Thus the MPFD could be a better candidate than design point for the mean of the sampling density function. This strategy will be investigated in the future work. To further examine their applicability, the SLP- and BND-based methods developed herein must be applied to practical cases for estimating lag-1 two-state transition probabilities and resultant resilience, particularly those with many random variables. Furthermore, the proposed SLP- and BND-based methods shall be extended for solving multi-lag problems, such as lag-2 resilience. 77 Chapter 3 develops an approach for seasonal evaluation of Pis in a water supply system with two or more rivers as the main water source. Such a system is rather common in practice. The approach should be applied to real cases in the future. 5.4 LIMITATIONS As approximating approaches, the SLP- and BND-based methods developed herein have limitations commonly faced by all other approximating techniques in reliability analysis. In particular, FORM-related drawbacks may affect the applicability of the two proposed methods. The linearization of highly irregular performance functions may lead to considerable errors. They may not readily apply to cases where an explicit performance function is not available. Although they may facilitate the estimate of resilience in cases where an explicit performance function is available, the methods do not apply when the performance function of a system is expressed as an implicit function. An implicit performance function is common in hydrotechnical systems, particularly in those where storage effects are important. For example, estimates of the optimal operation policy of a reservoir system for a given duration may require the analysis of a dynamic optimization program, which would consider accumulative effects of various benefits and costs. Thus Pis related to the operation policy will be characterized by an implicit performance function. For more information related to evaluating Pis for reservoirs, please see Hashimoto et al.(1982), Moy et al. (1986), and McMahon et al. (2006). 78 REFERENCES 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. McMahon, T. A., Adeloye, A. J., and Zhou, S.-L. (2006). "Understanding Performance Measures of Reservoirs." Journal of Hydrology, 324, 359-182. Moy, W. S., Cohen, 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. 79 APPENDIX A N U M E R I C A L APPROXIMATIONS O F DESIGN POINTS IN RELIABIL ITY ANALYSIS U N D E R P A R A M E T R I C C H A N G E S A version of this chapter has been accepted for publication in ASCE Journal of Engineering Mechanics. The paper is titled as Numerical Approximations of Design Points in Reliability Analysis under Parametric Changes by Y. Li and B. J. Lence. 80 PREFACE Traditional approaches for repeatedly updating reliability estimates, as needed in reliability-based optimal designs or real-time system control, require the iterative application of a reliability method. This paper explores a new strategy for repeatedly estimating reliability under frequent parameter variations. The central idea is to update the design point in the parameter domain, rather than in the traditional random variable domain, by evaluating several parametric sensitivity measures which are systems of nonlinear first-order ordinary differential equations relating the design point to parameter changes. Four numerical algorithms for evaluating the sensitivity measures are developed using the Euler and the improved Euler algorithms. Two solution procedures are applied. One procedure solves for the updated design point directly, while the other solves for both the unit normal vector and the radius of the design point separately, and evaluates the product of these to determine the updated design point. The numerical techniques are thoroughly compared with the classical HL-RF algorithm in five numerical examples regarding efficiency and accuracy. It is found that they are efficient and robust under given conditions. 81 A . l INTRODUCTION The reliability of civil engineering systems often must be evaluated under changing parameters due to changing conditions. (Parametric changes considered in this paper indicate any deterministic increment of the value of a set of parameters in a system, unlike a perturbation that normally represents an infinitesimal increment.) A changing parameter could be a distributional parameter of a random variable, i.e. mean or standard deviation under non-stationary conditions, or a deterministic parameter of a system's performance function. Deterministic parameters typically include decision variables to be determined by the designer or the decision maker, varying system properties such as changing environmental and anthological loads or aging material strength, as well as some factors for evaluating special risk-based indicators such as the resistance levels used In evaluating vulnerability, which indicates the likely magnitudes of a system failure in a naturally repairable system and may be estimated by evaluating the system's reliability under several different resistance levels (Hashimoto et al. 1982; Maier et al. 2001). In some cases, reliability must be repeatedly evaluated under very frequent changes of the system parameters. Typical examples are uncertainty-based optimal design and real-time system control problems. In the former, risk-based indicators of a system, i.e., reliability, may be incorporated in the constraints or in the objective function of a decision model that optimizes decision variables. Reliability must be frequently estimated given changing decision variables, whether a conventional optimization algorithm or a modern heuristic optimization technique such as a genetic algorithm is used to solve the model (Chen et al. 1999; Vasquez et al. 2000; Tolson et al. 2004). When vulnerability is considered, the reliability is also frequently estimated under different resistance levels. Furthermore, the repeated evaluation of reliability is indispensable in optimization models for time-variant risk problems in that the indicators may have to be evaluated over time considering changing distributional parameters due to non-stationarities. In the latter, real-time equipment control may require the repeated estimation of system reliability given frequently varying system parameter inputs. The estimates of system reliability may be useful in either indicating whether the equipment is approaching a malfunctioning state or in evaluating whether the frequently changing parameter values are indeed drawn from the same parent distribution. In these scenarios, of course, efficient updating of reliability estimates is desired. Practical estimating methods for reliability include simulations, approximations, and combinations of these (Ang and Tang 1984; Madsen et al. 1986; Sundararajan 1995; Tung et al. 2006). In particular, design point-based approximations, such as the first-order reliability method 82 (FORM) and the second-order reliability method (SORM), are widely used in practice because of their relative computational efficiency, an advantage over simulations, and their adaptability to complex problems, an advantage over the analytical method. An extremely important task in the approximations is the evaluation of the design point. The task is a nonlinear constrained optimization problem, mainly solved by optimization algorithms or iterative schemes (Liu and Der Kiureghian 1991). Optimization algorithms, i.e., the sequential quadratic programming (SQP) method or the gradient projection (GP) method, are usually robust, but require intensive evaluations of the performance function and its gradient, even when the gradient is analytically obtained. Here robustness indicates the ability to solve for the design point with sufficient accuracy (Liu and Der Kiureghian 1991). In addition, the complexity of the optimization algorithms affects their operability. One may not implement them properly without sufficient understanding of optimization programming, such as regarding search directions, step sizes, or merit functions. The iteration schemes, i.e., the Hasofer and Lind-Rackwitz and Fiessler (HL-RF) algorithm (Hasofer and Lind 1974; Rackwitz and Fiessler 1978), or the Ang-Tang method (Ang and Tang 1984), are relatively easier to operate, and more computationally efficient, requiring moderate evaluations of the performance function and its gradient. However, the convergence of these iteration schemes depends greatly on a proper starting point that is not easy to obtain in complex scenarios, thus they may not converge (Madsen et al. 1986; Liu and Der Kiureghian 1991; Tung 1996). To improve its robustness, Liu and Der Kiureghian (1991) and Zhang and Der Kiureghian (1997) have developed an improved HL-RF algorithm (iHL-RF) by introducing a line search scheme, thus increasing its complexity. If a design point-based approximation is employed in an environment where frequent estimates of reliability are required, reliability estimates can be updated by repeatedly applying the approximation. However, the repeated application may result in a high computational cost, and may also suffer from the risk of non-convergence.' Offering a potential alternative to the design point updating in environments where repeated estimates of reliability are needed due to frequent parametric changes, this paper explores in depth the feasibility of estimating the design point in the parameter domain, instead of in the traditional random variable domain, by numerically solving parametric sensitivity measures of the design point to parametric changes. Both domains are defined in the next section. This estimation approach was originally proposed by Bjerager and Krenk (1987; 1989). The use of this approach, enhanced by two algorithms and two solution procedures developed herein, is 83 explored using examples drawn from the reliability literature in terms of computational cost and effectiveness. The remaining sections of this paper are organized as follows. The design point and relative parametric sensitivity measures are first described. Then, numerical techniques are developed. Finally, the techniques are examined using several numerical examples. A.2 P A R A M E T R I C S E N S I T I V I T Y M E A S U R E S Reliability is defined based on a load-resistance interference characterized by a performance function g. This paper focuses on the component reliability problem where only one g exists. The probability of g > 0 is referred to as the reliability, and g - 0 defines a boundary called the limit-state surface. In FORM and SORM, original random variables are transformed into a standardized uncorrelated normal space (u) with zero means and unit standard deviations. The three approximations replace the original g with the first- or second-order Taylor Series of g expanded at the design point (u*). As a local optimum on the limit-state surface of g, u* is locally closest to the origin of the u-space. Here a basic assumption applies; that is, only one u* exists. Despite the existence of some scenarios with multiple design points (Gupta and Manohar 2004), the single-u* assumption has been proven effective in many applications. If g > 0 at the u origin, the reliability index /? equals the distance from the origin to the limit-state surface, or the norm of u*; otherwise, will be the negative of that distance. The two approaches approximate reliability using information regarding u*. In FORM, for example, reliability is approximated as 4>(y6), where <&() is the cumulative distribution function of the standard normal distribution. Also, u* provides important sensitivity information with respect to both u and the parameters. The unit normal vector to the limit state surface (or the negative of the normalized gradient vector of g) at u* a , indicates the sensitivity of f3 relative to the change of u*. (Note that u* = J3a, for u* is a local optimum.) Sensitivities of a , /?, and u* due to changes in parameters are determined by parametric sensitivity measures. The asymptotic approximation of the parametric sensitivity of f3 was originally developed by Hohenbichler and Rackwitz (1986). Then Bjerager and Krenk (1987; 1989) derived the following complete analytical expressions of the parametric sensitivity measures: du* dw (A.1) d w d/J |Vg dw (A.2a) 84 da _ 1 d w ~ / ? M ( _ * L - , » - A ^ i a + a T A / _ W t _ A ! L d w a Aa Vg d w Vg d w a Aa (A.2b) where w = the vector of desired parameters, i.e. deterministic parameters in g or distributional parameters of the original random variable; and Vg = the gradient vector of g relative to u at u* also known as the Jacobian. A is determined by A = ( I + ^ T H ^ ' <A-3> where I = an identity matrix; and H - the matrix of second-order partial derivatives of g relative to u at u*, also known as the Hessian. In fact, (A.2a) is the asymptotic approximation determined by Hohenbichler and Rackwitz (1986). These analytical expressions are derived based on another two basic assumptions in addition to one described above; that is, 1) g is twice differentiable relative to u, and 2) both g and Vg are differentiable relative to w . If a reliability-based optimization employs FORM, fj can serve as a substitute indicator for reliability, and thus (A.2a) is often used in such an optimization (Enevoldsen 1994; Pu et al. 1997). As a special realization of u, u* is a function of w , g, and the distribution of the original random variables. Conventional optimization algorithms or iterative schemes estimate u* in the u-space by searching among different u realizations, with fixed w . When w changes, these methods repeat the searching process with a changed w . In this case the u-space is defined as the random variable domain. Since (A.l) and (A.2) are exact derivatives of w , it is likely that, under certain conditions, u* can be updated given changing w by solving for these sensitivity measures, without the need to search in the random variable domain. In other words, u* is estimated in the w-space that is defined as the parameter domain, as shown in the next section. A.3 NUMERICAL APPROXIMATION ALGORITHMS Each of the equations in (A.l) and (A.2) is a system of nonlinear first-order ordinary differential equations (ODEs). It is nearly impossible to integrate these three interrelated measures in practice due to their nonlinearity. Numerical techniques for solving them are herein proposed. The techniques apply to a component reliability problem with a set of standardized random variables, u (u\, u2,..., um)r, and a set of parameters, w (w,, w 2 , . . . , w„) T, which satisfies three conditions; that is, 1) g is twice differentiable relative to u, 2) only one u* exists, and 3) both g and Vg are differentiable relative to w . 85 A.3.1 Numerical schemes A critical strategy for developing the numerical schemes is to reduce the numerical problem to an initial-value problem (Figure A . l ) . A reference design point, u0, is specified based on an initial parameter set, wo, by any existing optimization algorithms or iterative schemes. Given a new parameter set W i , the numerical schemes approximate the design point for wi (u]), by solving the ODEs with an initial value (u*Q). For efficiency purposes, low-order numerical techniques should be adopted wherever possible. This paper adopts the Euler and the improved Euler algorithms to solve the ODEs. To further enhance efficiency, for each parametric change, the step length vector, h (hi, h2,..., hn)T, is set to be W i - wr> Thus stability and convergence problems are eliminated. J Two solution approaches may be applicable for (A.l) and (A.2). After substituting (A.2a) and (A.2b) into (A.l), u* can be updated by simply solving for (A.l). This approach is called the joint solution procedure. Alternatively, (A.2a) and (A.2b) are solved separately, and then the product of the solutions is used to update u* given that u* = J3 a. This second approach is called the separate solution procedure. Here the separate solution procedures with the Euler and the improved Euler algorithms are denoted as S-Euler and S-iEuler, respectively, and the joint solution procedures with the Euler and the improved Euler algorithms as J-Euler and J-iEuler, respectively. J-Euler and J-iEuler solve the paired sets of equations Figure A . l . Configuration of initial-value problem in a bivariate u-space u P = u, + h T du* (A.4a) dw 86 * * N T r d u * * d u * * u c =u 0 +0.5h [ - — ( w 0 , u 0 ) + - — ( w „ u p ) ] (A.4b) dw dw where u * = a predictor of u]; and u^, = a corrector for u * . In the paired sets, (A.4a) is the solution of J-Euler, while (A.4b) is the solution of J-iEuler. S-Euler and S-iEuler solve the following sets of equations . ^=/?0 + h T M ( W o ) / ? o ) ( A . 5 a ) dw a P = a 0 + h T - ^ - ( w 0 , a 0 ) (A.5b) d w u* = ppav (A.6) pc = / 5 0 + 0 . 5 h T [ ^ ( w 0 , ^ 0 ) + M ( w 1 , ^ ) ] (A.7a) dw dw a c =a0+0.5hT[4^-(w0,a0)+-^-(w„aP)] (A.7b) d w d w U c = / ? c a c (A.8) where p\ = the reliability index for u * ; fjP = a predictor of the reliability index for u*; ao = the unit normal vector at u j ; a p = a predictor of the unit normal vector at Uj ; pc = a corrector for PP; and Oc = a corrector for ap. Note that (A.5b) and (A.7b) must be normalized in real applications. In the system, (A.6) is the solution of S-Euler, while (A.8) is the solution of S-iEuler. Equations (A.5a) and (A.5b) have been recommended by Bjerager and Krenk (1987; 1989) and Santos et al. (1995) as a rapid approximation technique. As design point is a vector, the error of a numerical technique can be measured as the absolute value of the distance between u* and the solution of the technique. (Note that this error measure differs from the difference between two reliability indexes, and is commonly larger than that difference.) Since h = wi - wo, no sub-steps are applied in the numerical techniques. As a result, all numerical techniques presented above exhibit only local error, and no global errors are created. Ignoring round-off errors, elementary knowledge of numerical methods for solving ODEs indicates that the local error of the improved Euler algorithm is of third order, i.e. o(||h||3), while the Euler algorithm has a second-order local error, i.e. o(||h|| 2). On the other hand, due to the crude definition of step length, the third-order error may not be necessarily smaller than the second-order error; that is, in some cases an improved Euler algorithm may be less accurate than 87 an Euler algorithm. In the mean time, the errors of both procedures are of the same order; that is, both S-Euler and J-Euler have a second-order error, while both S-iEuler and J-iEuler have a third-order error. It is theoretically difficult to evaluate the most accurate procedure a priori. So the most accurate numerical technique is likely to be case specific. The crude definition of step length might lead to a large probable maximum error (PME) in a given case. The PME usually occurs at one of the corners of the w domain; that is, where the values of the w components reach one of their combined limits. In civil engineering systems, however, parametric values are limited, and consideration of systems experiencing widely varying changes of parametric values is not common. Decision variables are typically constrained by technical codes, by physical constraints, or by economical constraints. Reliability analyses of systems facing excessive parametric changes would likely indicate that the system is either too safe or very risky and therefore should be excluded from normal cost-safety considerations. Hence, a reasonably moderate range of possible values should be considered for each component. This will have an effect of controlling the PME. A.3.2 Computational procedures Computational procedures of the techniques are described using the working flowcharts in Figures A.2 and A.3, and computational costs are roughly estimated based on these. The cost budget in this paper merely considers the computational burden incurred while updating u* for a parametric change, reflecting the characteristic of an environment featured by frequent parametric changes. The computational cost is represented by the number of substantial terms evaluated during the updating process, and is referred to as the cost index (CI). The working flowchart for J-iEuler (Figure A.2) contains eight steps, divided into the preparation and the operation procedures. The preparation procedure provides basic information needed to initiate a reliability-updating task. While costly, the preparation procedure needs to be undertaken only once, just prior to the updating process. Therefore the cost of this procedure is excluded from the cost budget. The operation procedure is undertaken for each parametric change, and is considered in the cost budget. Step 5 solves J-Euler. The computational cost of an Euler algorithm is negligible (Powers 1987), particularly considering that Step 5 utilizes known information from the preparation procedure. Similarly, Step 7 requires little computational effort. The first two terms in Step 6 are easily obtained from u P , so their evaluation is inexpensive. The most computationally expensive effort is undertaken to evaluate the last five terms at Step 6. When all derivatives involved at Step 6 are analytically computed, the CI of J-iEuler is therefore 88 five. If these derivatives are numerically computed using a central difference scheme, the number of evaluations of g during an updating process can be used to determine the CI. Specifically, Vg and H require evaluating g for 2m2 + 1 times, and dg/dw and 3Vg/3w for Amn + 2n times. It is difficult to quantify the expense of the evaluation of A, which depends on its size and solution algorithm. In this case, the CI of the J-iEuler technique is hence 2m + Amn + 2n + 1 plus the equivalence of A. It is clear that J-Euler stops at Step 5, incurring a negligible cost. Preparation Operation Step 1: Define the range of w i Step 2: Specify w 0 Step 3: Find u 0 based on W Q Step 4: Evaluate fi, Ot, Vg, H, A , dg/dw, and dVg/dw using W Q and u„ Step 5: Compute up using (A.4a) Step 6: Evaluate fi, a, Vg, H, A , dg/dw, and dVg/dw using up and W ) I Step 7: Compute u*c using (A.4b) Figure A.2. The working flowchart of J-Euler and J-iEuler The working flowchart of S-Euler and S-iEuler has the same preparation procedure as that of J-iEuler. The operation procedure of the S-iEuler technique is illustrated in Figure A.3. The substantial effort in this technique also is undertaken at Step 6, where five important terms are evaluated. Consequently, the CI of S-iEuler is the same as that of J-iEuler whether the derivatives involved are analytically or numerically computed. S-Euler simply solves Equation (A.6) using information provided at Step 5, and a negligible cost is incurred. For comparison sake, the computational cost of the classical HL-RF algorithm is roughly estimated as follows. HL-RF might be the fastest algorithm for searching for design point in practice. Considering the substantial effort required for evaluating Vg and g, the CI per iteration of HL-RF would be two when Vg is analytically computed, or 2m +1 when Vg is numerically computed. 89 Step 5.1: Compute BP using (A.5a) 3 Step 5.2: Compute a P using (A.5b) Step 6: Evaluate Vg, H , A, dg/dvt and dVg/dw using W! and u p Step 7.1: Compute fic using (A.7a) Step 7.2: Compute Oc using (A.7b) Step 8: Compound u *c using (A.8) Figure A . 3 . The operation procedure of S-iEuler If all required derivatives are analytically computed, J-iEuler or S-iEuler will be approximately equivalent to three iterations of HL-RF in terms of CI, while J-Euler or S-Euler causes a negligible cost. Hence the parameter domain-based numerical techniques may be compared with HL-RF in terms of efficiency through comparing their errors, as illustrated in Figure A . 4 . Assuming that they all start from u * , for example, if the error of HL-RF at the first iteration is smaller than J-Euler or S-Euler, HL-RF will be relatively faster than these algorithms. Likewise, if the error of HL-RF at the third iteration is larger than J-iEuler or S-iEuler, HL-RF will be relatively slower than these algorithms. Figure A . 4 . Illustration of efficiency comparison between HL-RF and numerical techniques 90 Sometimes g is expressed in a complicated form, for example as a multivariate and highly nonlinear function, or just implicitly, so all derivatives required in a reliability analysis might have to be numerically computed, for instance by a finite difference scheme, leading to expensive workloads in implementing the numerical techniques. An adaptation for this case is to ignore second-order derivatives, H and dVg/dw, as suggested by Bjerager and Krenk (1987; 1989). The approach could cause considerable errors, so in a given case one must investigate its justification before adopting it. If the adaptation is not justified, J-iEuler and S-iEuler will generally lead to a high computational cost, indicated by their CI, and thus are not recommended. In a given case, if an Euler algorithm is sufficiently accurate, as identified by the tolerance requirements for Karush-Kuhn-Tucker (KKT) conditions or for probabilistic contents, however, it may be worth applying the Euler technique, because of its quite small computational cost in the operation procedure, despite numerical efforts spent at Step 4 as in Figure A.2. Note that Step 4 is undertaken only once for the entire reliability-updating task, while HL-RF needs to be run for each parametric change during the task. Although tedious and burdensome, a finite difference scheme is normally reliable and feasible (Haldar and Mahadevan 1995). The efficiency comparison implies that a combination of the numerical techniques and HL-RF may facilitate updating design point. If more efficient than HL-RF for limit combinations of parameter values, a numerical algorithm can be used to update reliability under quickly changing parameters. Upon each parametric change during the task, this technique is run first to produce an initial estimate. If it meets the design point conditions, namely KKT conditions, this initial estimate just is the exact new design point; otherwise one can run HL-RF starting from the initial estimate to search for the exact solution. The combined application accelerates HL-RF by saving iterations, as shown in two numerical examples in the next section. A . 4 N U M E R I C A L E X A M P L E S To demonstrate the feasibility of the parameter domain-based numerical techniques, five numerical examples are examined. The examples are designed to represent a screening process that identifies a robust and efficient method under frequent parametric changes. All of the numerical techniques and HL-RF are applied to the limit combinations of parameters, thus producing their PMEs. The numerical techniques are compared with HL-RF in terms of efficiency. Specifically, the errors at the first and third iterations of HL-RF are compared with the errors of the numerical techniques. In each example, a comprehensively effective method is ( chosen as the tool for updating design point under parametric changes. 91 In al l examples, Wo is set to be the median o f w. T o compare fairly, H L - R F searches for u] starting from u„, without finding a better starting point. A l l required derivatives are analytically obtained. The convergence criterion of H L - R F adopts an absolute tolerance of 10"8. A.4.1 Example 1: A plane frame structure The combined failure mechanism of a plane frame structure shown in Figure A . 5 is characterized by the fo l lowing g: g = Zl+2z3 + 2 z 4 + z5 - FH Lh - FV Lw (A.9) where zi,---, zs = the plastic moment capacities ( k N m); FH = the horizontal load (kN); FV = the vertical load (kN); L/, = the height o f the structure (m); and Lw = the half-span length o f the structure (m). This example is taken from Example 5.6 in Madsen et al . (1986). Suppose that the changing parameters are Lh and Lw, each o f which lies in the range [4.5, 5.5]. The other terms are mutually statistically independent lognormal random variables, whose statistical properties are summarized in Table A . l . Figure A . 5 . A plane frame structure The H L - R F algorithm produces u 0 starting from the origin of the u-space (J3n = 2.8825). Four combinations of parametric changes are evaluated. H L - R F converges for each combination starting from u0, and the resultant solution is considered as u,. The results, as shown in Table A . 2 , provide 1) the errors of the numerical techniques, 2) the errors of H L - R F at the first and third iterations, 3) the total number o f iterations taken by H L - R F to search for u*, and 4) the reliabili ty index for u*, or j3\. 92 Table A. 1. Statistical properties of random variables in numerical examples Example Random Mean Standard variables deviation 1 Zl , . . . , z5 134.9 kNm 13.49 kNm FH 50 kN 15 kN Fv 40 kN 12kN 2 T, 360 MPa 36 MPa T 1 c 40 MPa 6 MPa D 0.3 m 0.015 m 3 250 kNm 75 kNm Xl 125 kNm 37.5 kN m x-i 2500 kN 500 kN X4 40 MPa 4 MPa 4 V 16.09 km/d 4.827 km/d Kd 0.35 /d 0.10/d Ka 0.7 Id 0.2 Id 18 mg/L 5 mg/L D0 1.0 mg/L 0.3 mg/L 5 X\,..., X4 120 12 xs 50 15 x6 40 12 Table A.2. Computational results for Example 1: A plane frame structure Parameters Errors of proposed numerical techniques HL-RF Error at Error at Total No. J-Euler J-iEuler S-Euler S-iEuler the 1 st the 3rd number of (m) (m) fii iteration iteration iterations 1 5.5 5.5 2.4477 0.0351 0.0024 0.0240 0.0008 0.0993 0.0168 10 2 4.5 5.5 2.9715 0.1581 0.0089 0.1293 0.0138 0.3442 0.1212 16 3 4.5 4.5 3.3588 0.0443 0.0027 0.0307 0.0013 0.1788 0.0461 13 4 5.5 4.5 2.7603 0.1158 0.0027 0.0867 0.0041 0.1365 0.0201 9 All numerical techniques appear to be more efficient than HL-RF. S-Euler is superior to J-Euler, while J-iEuler and S-iEuler tie in general. J-iEuler and S-iEuler obviously improve J-Euler and S-Euler, respectively. As a promising reliability-updating method in this example, J-iEuler can be safely used under an error level of 0.009. 93 To illustrate its accelerating effect, the combined application of the numerical techniques and HL-RF is implemented on the four combinations. Assume that S-Euler is to be used in the combined application, and its solutions have not met KKT conditions. (Note that meeting KKT conditions depends on user-specified tolerance criteria in reality.) Starting from the S-Euler solutions, HL-RF iterates merely 7, 13,9, and 8 times to converge for parametric combinations 1 to 4, respectively. The net savings of iteration are 3, 3, 4, and 1, respectively. In conjunction with a parameter domain-based numerical method, the efficiency of HL-RF is thus improved. A.4.2 Example 2: A reinforced concrete beam This example is developed based on Example 4.3 in Madsen et al. (1986). The bending moment safety of a reinforced concrete beam is evaluated by the following g: KT2 A 2 8 = DTsAs f-^-Mb (A.10) where Mb - the sectional bending moment (MNm); As = the area of reinforcement (m ); Ts = the yield stress of the reinforcement (MPa); Tc- the maximum compressive concrete strength (MPa); Sb = the cross-sectional width (m); D = the effective depth of the reinforcement (m); and K = a factor dependent on the stress-strain relation of concrete. Suppose that Sb and As change in [0.1, 0.15] and [200x10"6, 250xl0"6], respectively. Assume Mb and K to be deterministic, 0.016 and 0.5, respectively. The other terms are mutually independent normal variables whose statistical properties are provided in Table A.l . The HL-RF algorithm produces uj starting from the u origin (j5o = 3.0894). As shown in Table A.3, compared with HL-RF, S-iEuler is obviously more efficient except for one combination, and the other numerical techniques are also competitive. In this example, S-iEuler can be an effective reliability-updating tool, and its error in estimating reliability will not exceed 0.000035 in FORM. It is noted that both S-Euler and J-Euler cause an error of about 0.003 in the reliability estimate by FORM. If this error of reliability estimate is acceptable in this example, both may be used as an efficient reliability-updating tool with a small cost. 94 Table A.3. Computational results for Example 2: A reinforced concrete beam Parameters Errors of proposed numerical techniques HL-RF Error at Error at Total As Sb No. J-Euler J-iEuler S-Euler S-iEuler the 1st the 3rd number of (10'6m2) (m) p\ iteration iteration iterations 1 250 0.15 3.7808 0.0774 0.0013 0.0743 0.0004 0.1220 0.0021 6 2 200 0.15 2.2966 0.0873 0.0026 0.0850 0.0003 0.0749 0.0005 5 3 200 0.1 2.2331 0.0962 0.0042 0.0916 0.0010 0.0825 0.0005 5 4 250 0.1 3.7251 0.0740 0.0013 0.0731 0.0005 0.1107 0.0021 6 A.4.3 Example 3: A rectangular column The following g is used for designing the rectangular column in Example 4.3 of Royset et al. (2001): 8 = 1 - ^ - - - ^ — i - ^ - ) 2 (A.11) where Sh= the cross-sectional height (m); and x i , . . . , x 4 = statistically independent lognormal variables, see Table A.l for their properties. Assume that Sb and Sh range in [0.3, 0.4] and [0.5, 0.6], respectively. The results are shown in Table A.4. The HL-RF algorithm produces u* starting from the u origin (fio - 3.0506). The numerical algorithms basically outperform HL-RF. In this example, S-iEuler may be chosen as a reliability-updating tool, and its error in estimating reliability will not exceed 0.001 using FORM. Table A.4. Computational results for Example 3: A rectangular column Parameters Errors of proposed numerical techniques HL-RF No. sb (m) s,, (m) Px J-Euler J-iEuler S-Euler S-iEuler Error at the 1st iteration Error at the 3rd iteration Total number of iterations 1 0.4 0.6 4.6013 0.4981 0.0593 0.3722 0.0881 0.6804 0.2794 26 2 0.3 0.6 2.6046 0.2163 0.0188 0.1683 0.0095 0.1760 0.0252 10 3 0.3 0.5 1.2562 0.2373 0.1390 0.1296 0.0052 0.3379 0.0062 7 4 0.4 0.5 3.2014 0.1835 0.0068 0.1506 0.0091 0.2391 0.0648 13 Starting from the solution produced by S-Euler, HL-RF iterates 24, 9, 6, and 7 times to converge for combinations 1 to 4, respectively; the net savings in terms of iterations are 2, 1, 1, 95 and 6, respectively. Starting from the solution produced by S-iEuler, HL-RF iterates 18, 6, 3, and 7 times to converge for combinations 1 to 4, respectively; the net savings are 5, 1, 1, and 3, respectively. A.4.4 Example 4: Riverine dissolved oxygen response In a riverine dissolved oxygen (DO) response system governed by the Streeter-Phelps model (Streeter and Phelps 1925), the reliability that, at a given location downstream of a point source discharge of carbonaceous biochemical oxygen demand (CBOD) waste, the DO deficit remains below an allowable deficit level is characterized by the following g: g=Da £ t * 2 _ ( g - * , . v / v - e - W ^ - D ^ - W (A-12) Ka~Kd where Da = the allowable DO deficit (mg/L); y'= the distance from the source (km); v = the average flow velocity (km/d); K~d = the deoxygenation coefficient (1/d); Ka = the reaeration coefficient (1/d); Lo = the initial CBOD concentration immediately downstream of the point source (mg/L); and Do = the initial DO deficit immediately downstream of the point source (mg/L). Assume that v, Kj, Ka, Lo, and DQ are mutually statistically independent lognormal variables, whose statistical properties are provided in Table A.l. Suppose that Da and y are in the range [2, 8] and [8.045, 16.09], respectively. In this reliability problem with a specified y, FORM has been proven effective for a wide range of values of j5 (Tung 1990; Melching and Anmangandla 1992). The HL-RF algorithm produces u „ starting from the u origin [fio - 0.8965), and converges for each parametric change starting from u „ , thus producing u\. The results (Table A.5) show that all numerical schemes except J-iEuler outperform HL-RF for two combinations. In this example, no single method is apparently better than the others. Table A.5. Computational results for Example 4: Riverine dissolved oxygen response Parameters Errors of proposed numerical techniques HL-RF Da (mg/L) y (km) Error at Error at Total No. J-Euler J-iEuler S-Euler S-iEuler the 1st iteration the 3rd iteration number of iterations 1 8 16.09 1.9232 0.3092 0.0328 0.2675 0.0051 0.2217 ' 0.0070 9 2 2 16.09 -2.0984 0.9744 0.6054 0.8979 0.2004 1.1649 0.0388 10 3 2 8.045 -1.2735 0.9415 0.8719 0.8891 0.1783 0.7245 0.0016 6 4 8 8.045 2.7440 0.3060 0.0387 0.2687 0.0204 0.8882 0.0389 9 96 A.4.5 Example 5: Impacts of noise The performance function in Example 2 of Liu and Der Kiureghian (1991) tests the impacts of noise on system performance. It contains a high-frequency artificial noise term: 6 g = x,+2x2+2x3+x4-5x5-5x6-l-/J^sin(100x/) (A. 13) /=i where x\,..., X(, = mutually independent lognormal variables, whose statistical properties are shown in Table A.l; and B = a constant in the range [0, 0.001]. In this highly nonlinear problem, HL-RF exhibits unstable behaviors. It does not converge after 100 iterations for producing u0 starting from the u origin. Thus an optimization algorithm is used to produce u j (fio = 2.3482). The results are shown in Table A.6. When 6 = 0, the HL-RF algorithm converges starting from u j , and the resultant solution is considered as u j . When B = 0.001, however, HL-RF does not converge starting from u j after 100 iterations, implying a weak robustness. Several optimization algorithms and iterative schemes, GP, SQP, and iHL-RF, converge starting from the u origin, and their solution is considered to be the actual solution of u , . The Cis of GP, SQP, and iHL-RF are 236, 275, and 242, respectively, evaluated based on Liu and Der Kiureghian (1991). Although more accurate, these algorithms and schemes are much more computationally expensive than the numerical techniques. In this example, both S-Euler and J-Euler are a good choice as the reliability-updating tool with a PME of 0.01, as both are relatively more robust than HL-RF, and more efficient than S-iEuler, J-iEuler, GP, SQP, and iHL-RF. Under this PME, the error of reliability estimate caused by .S-Euler or J-Euler will not exceed 0.0003 in FORM. Despite almost identical reliability indexes, this example comprehensively demonstrates the robustness and efficiency of the parameter domain-based numerical techniques. Table A.6. Computational results for Example 5: Impacts of noise Errors of proposed numerical techniques HL-RF Error at Error at Total B p\ J-Euler J-iEuler S-Euler S-iEuler the 1st the 3rd number of iteration iteration iterations 0 0.001 2.3482 2.3482 0.0058 0.0103 0.0011 0.0101 0.0052 0.0105 0.0039 0.0104 0.0002 0 not converg 3 ent 97 A .5 C O N C L U S I O N S The practice of civil engineering reliability analysis over the past three decades has shown that no single method is perfect in searching for u*, in terms of efficiency and robustness. This is why professionals have developed so many methods. Given the need to frequently estimate reliability, as encountered in uncertainty-based optimal designs and real-time system control problems, it will be meaningful to identify a robust and efficient reliability-updating method for a specific case. Numerical techniques based on Bjerager and Krenk (1987; 1989) are thus investigated in depth here as alternatives, whereby u* is estimated within the parameter domain using parametric sensitivity measures that comprise several systems of complicated nonlinear first-order ODEs. The parameter domain-based numerical techniques reflect the characteristics of an environment where parameters are frequently changing. To solve these ODEs practically, four numerical schemes have been developed based on the Euler and the improved Euler algorithms, as well as two solution procedures. Especially, S-Euler and J-Euler require a negligible computational cost. The feasibility and behaviors of the numerical methods are demonstrated in numerical examples where required derivatives are analytically obtained. In most cases, the separate procedure is more accurate than the joint procedure, and the improved Euler algorithms improve on the Euler algorithms in this reliability-updating problem. One can compare the efficiencies of the numerical techniques and HL-RF through comparing their errors related to the limit combinations of parameters. It is observed that the numerical techniques seem more efficient than HL-RF in some cases. A combined application of the numerical techniques and HL-RF has been proposed, and the reliability updating may accelerate in this combined application. It must be noted that the applicability of any numerical technique presented in this paper is case-specific, depending on the actual situation. Generally speaking, in screening a robust and efficient reliability-updating tool adaptive for frequent parametric changes, the parameter domain-based method should be seriously considered for practical applications when 1) iterative schemes, especially the HL-RF algorithm, are not convergent, or 2) iterative schemes are convergent, but less efficient than the numerical methods. Nevertheless, the parameter domain-based numerical methods involve evaluating several first- and second-order derivatives, which are often numerically computed in complex situations. In this case, S-Euler and J-Euler are still promising if proven sufficiently accurate. The CI of both S-iEuler and J-iEuler significantly increases as m or n increases, limiting their applicability. 98 (While ignoring second-order derivatives largely decreases computational efforts, the justification of such ignorance is case-based.) Thereafter it is recommended that S-iEuler or J-iEuler be considered when all or a part of required derivatives can be analytically obtained. For example, in the reliability-based optimization design of a riverine point-source water quality management program (Vasquez et al. 2000), if the steady-state DO concentration is estimated by the Streeter-Phelps model or its modifications (Camp 1963; Dobbins 1964), the DO concentration at any location in a river is a linear function of the decision variables of the problem, namely the treatment efficiencies of upstream wastewater treatment plants. As a result, g is a linear function of w, basically avoiding the numerical work required to estimate dg/dw and dVg/dw, and S-iEuler or J-iEuler might thus be worth testing. 99 References Ang, A. H. S., and Tang, W. H. (1984). Probability Concepts in Engineering Planning and Design, Vol. II: Decision, Risk, and Reliability, John Wiley, New York. Bjerager, P., and Krenk, S. (1987). "Sensitivity Measures in Structural Reliability Analysis." The First JFIP Working Conference on Reliability and Optimization of Structural Systems, P. Thoft-Christensen, ed., Springer-Verlag, 459-470. Bjerager, P., and Krenk, S. (1989). "Parametric Sensitivity in First Order Reliability Theory." ASCE Journal of Engineering Mechanics, 115(7), 1577-1582. Camp, T. R. (1963). Water and Its Impurities, Rheinhold Publishing, New York, N.Y. Chen, J., Wang, F., Sun, H., and Zhang, C. (1999). "A Method of Optimum Design Based on Reliability for Antenna Structures." Structural Engineering and Mechanics, 8(4), 401-410. Dobbins, W. E. (1964). "BOD and Oxygen Relationships in Streams." ASCE Journal of Sanitary Engineering Division, 90(3), 53-78. Enevoldsen, I. (1994). "Sensitivity Analysis of Reliability-Based Optimal Solution." ASCE Journal of Engineering Mechanics, 120(1), 198-205. Gupta, S., and Manohar, C. S. (2004). "An Improved Response Surface Method for the Determination of Failure Probability and Importance Measures." Structural Safety, 26, 123-139. Haldar, A., and Mahadevan, S. (1995). "First-Order and Second-Order Reliability Methods." Probabilistic Structural Mechanics Handbook, C. Sundararajan, ed., Chapman & Hall, New York, 27-52. 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. C. (1974). "Exact and Invariant Second-Moment Code Format." ASCE Journal of the Engineering Mechanics Division, 100(1), 111-121. Hohenbichler, M., and Rackwitz, R. (1986). "Sensitivity and Importance Measures in Structural Reliability " Civil Engineering Systems, 3(4), 203-209. Liu, P.-L., and Der Kiureghian, A. (1991). "Optimization Algorithms for Structural Reliability." Structural Safety, 9, 161-177. Madsen, H. O., Krenk, S., and Link, N. C. (1986). Methods of Structural Safety, Prentice-Hall, Englewood Cliffs, N.J. Maier, H. R., Lence, B. J., Tolson, B. A., and Foschi, R. O. (2001). "First-Order Reliability Method for Estimating Reliability, Vulnerability and Resiliency." Water Resources Research, 37(3), 779-790. Melching, C. S., and Anmangandla, S. (1992). "Improved First-Order Uncertainty Method For Water-Quality Modeling." ASCE Journal of Environmental Engineering, 118(5), 791-805. Powers, D. L. (1987). Elementary Differential Equations, PWS Publishers, Boston. Pu, Y., Das, P. K., and Faulkner, D. (1997). "A Strategy for Reliability-Based Optimization." Engineering Structures, 19(3), 276-282. Rackwitz, R., and Fiessler, B. (1978). "Structural Reliability under Combined Random Load Sequences." Computers & Structures, 9, 489-494. Royset, J. O., Der Kiureghian, A., and Polak, E. (2001). "Reliability-Based Optimal Design of Series Structural Systems " ASCE Journal of Engineering Mechanics, 127(6), 607-614. 100 Santos, J. L. T., Siemaszko, A., Gollwitzer, S., and Rackwitz, R. (1995). "Continuum Sensitivity Method for Reliability-Based Structural Design and Optimization." Mechanics of Structures and Machines, 23(4), 497-520. Streeter, H. W., and Phelps, E. B. (1925). "A Study of The Pollution and Purification of The Ohio River, III, Factors Concerned in The Phenomena of Oxidation and Reaeration." Public Health Bulletin 146, U.S. Public Health Service, Washington, D.C. Sundararajan, C. (1995). "Probabilistic Structural Mechanics Handbook." Chapman & Hall, New York. Tolson, B. A., Maier, H. R., Simpson, A. R., and Lence, B. J. (2004). "Genetic Algorithms for Reliability-Based Optimization of Water Distribution Systems." ASCE Journal of Water Resources Planning and Management, 130(1), 63-72. Tung, Y.-K. (1990). "Evaluating the Probability of Violating Dissolved Oxygen Standard." Ecological Modelling, 51(3-4), 193-204. Tung, Y. K. (1996). "Uncertainty and Reliability Analysis." Water Resources Handbook, L. W. Mays, ed., McGraw-Hill, 7.1-7.65. Tung, Y. K., Yen, B. C , and Melching, C. S. (2006). Hydrosystems Engineering Reliability Assessment and Risk Analysis, McGraw-Hill, New York. Vasquez, J. A., Maier, H. R., Lence, B. J., and Tolson, B. A. (2000). "Achieving Water Quality System Reliability Using Genetic Algorithms." ASCE Journal of Environmental Engineering, 126(10), 954-962. Zhang, Y., and Der Kiureghian, A. (1997). "Finite Element Reliability Methods for Inelastic Structures " UCB/SEMM-97/05, Department of Civil of Environmental Engineering, University of California, Berkeley, CA. 101
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Estimating resilience for hydrotechnical systems
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Estimating resilience for hydrotechnical systems Li, Yi 2007
pdf
Page Metadata
Item Metadata
Title | Estimating resilience for hydrotechnical systems |
Creator |
Li, Yi |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | With the cycle of natural conditions (hydrological or meteorological), a repairable hydrotechnical system may have a recovery capacity after a system failure. Resilience is a probabilistic performance indicator that measures such a recovery capacity. It is also one of performance indicators that can be used to measure the sustainability of a water resources project. In nature, resilience is a conditional probability that, given a system failure at an earlier time moment, the system can resume functioning at a later time moment. The conditional probability is a close derivative of two-state transition probabilities in a stochastic process. The current estimating techniques for the two-state lag-I transition probabilities may require a high computational cost, or a complex implementation procedure. This thesis explores practical approximation methods for estimating the two-state lag-1 transition probabilities in discrete processes, and the resultant resilience. Two methods are proposed herein. By developing a concise and straightforward approximation for the lag-1 autocorrelation coefficient of system performance functions at two consecutive moments on the basis of the first-order reliability method, one method improves the conventional approaches that estimate the transition probabilities based on a bivanate normal distribution of the system performance functions. The other employs a linear stochastic prediction of the performance function based on the mean point of failure or safety domain. The mean points are one of the fundamental properties of reliability problem. The thesis mathematically defines the mean points, and develops a practical approach for approximating them based on the first-order reliability method. The mean points-based method improves the conventional approaches that estimate the transition probabilities based on a linear stochastic prediction of performance function. The two methods developed herein are demonstrated in a typical river hydrology example for estimating seasonal lag-1 resilience. The proposed approximation methods for estimating the transition probabilities in discrete processes can be also applicable to other civil engineering problems, such as marine/offshore steel structure fatigue analysis, river navigability study, or concrete reinforcement corrosion prediction. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-16 |
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. |
IsShownAt | 10.14288/1.0063241 |
URI | http://hdl.handle.net/2429/31384 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-318140.pdf [ 5.51MB ]
- Metadata
- JSON: 831-1.0063241.json
- JSON-LD: 831-1.0063241-ld.json
- RDF/XML (Pretty): 831-1.0063241-rdf.xml
- RDF/JSON: 831-1.0063241-rdf.json
- Turtle: 831-1.0063241-turtle.txt
- N-Triples: 831-1.0063241-rdf-ntriples.txt
- Original Record: 831-1.0063241-source.json
- Full Text
- 831-1.0063241-fulltext.txt
- Citation
- 831-1.0063241.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-0063241/manifest