On reservoir reliability

Journal of Hydrology, 47 (1980) 269--296

269

Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands

[3]

ON R E S E R V O I R RELIABILITY

G.G.S. PEGRAM

Department of Civil Engineering, University of Natal, Durban 4001 (South Africa) (Received October 10, 1979; accepted for publication October 19, 1979)

ABSTRACT Pegram, G.G.S., 1980. On reservoir reliability. J. Hydrol., 47: 269--296. Finite-difference equations, integral equations and simulation are employed to ascertain mean first passage times (in particular the mean recurrence time of emptiness) of a non-seasonal reservoir fed by various input distributions. These include the normal, lognormal and discrete -+1 inputs, both independent and serially correlated. Included are descriptions of each method and some watchpoints. The results are accurate to within three significant figures or better; it was found that for large capacity, skewness, correlation and mean net input, great difficulty was experienced in obtaining accurate results. The concept of the standard reservoir (capacity expressed as a ratio to the standard deviation of the net input) facilitates comparisons between remarkably disparate marginal distributions of input. Results on discrete reservoirs yield respectable approximations to the continuous state--space reservoir.

INTRODUCTION

The reliability of a reservoir is the degree of dependability of its water supply for domestic, industrial or agricultural use. C o m m o n l y used measures of reliability o f reservoirs are the "probability of em pt i ness" or alternatively the "probability t hat de m and exceeds supply". Related measures of reliability are " t h e time to fill" or the " t i m e to e m p t y " . More accurately, these should be defined as first passage times, for then we can c o m p u t e measures such as " t h e mean first passage time from e m p t y (or full) to full (or e m p t y ) " . These numerical quantities are then used in the design of reservoir supply systems. The problem is to c o m p u t e them, and we look at three m e t h o d s in this paper: difference equations, integral equations and simulation. The problem of determining the probability distribution of storage can be posed as one of finding the solution to an integral equation in discrete time (Moran, 1959, p. 40). The most desirable solutions are analytical, but these are difficult to come by. The reason t hat t hey are desirable is because the quantities o f interest can be c o m p u t e d to an arbitrary degree of accuracy, when th ey exist in closed form. (Of course not all analytical results are

0022-1694/80/0000--0000/$02.25

O 1980 Elsevier Scientific Publishing Company

270 simple functions, and some, like the integral of the normal probability density function, must be computed from approximating series solutions.) Essentially, then, one can talk about "accurate" solutions, accurate to within six significant figures say, and " a p p r o x i m a t e " solutions where accuracy is measured to within about three significant figures, even when large amounts of computing effort have been used to acquire the figures. (Of course "app r o x i m a t e " solutions tend to become more "accurate" as one acquires greater computing power.) The methods used herein fall into three groups: (1) those based on the difference equation derived from the integral equation; (2) methods of solving the integral equation directly; and (3) simulation. Numerical approximations to the differential equations defining reservoir behaviour can be obtained in two ways. The one way is to define a system of linear difference equations in terms of probabilities. Alternatively, one may be able to solve a similar discrete state--space problem (which approximates the original problem) exactly. Both these approaches will be exploited in this paper. A second approach to solving the storage problem is by solving the integral equation directly. Analytical solutions are available in a few special cases, but useful approximations with known error bounds can be c o m p u t e d by well-documented numerical methods, especially in the independent input case. In this paper, we use integral equations to obtain solutions to corroborate some of the difference equation solutions. Falling into the class of approximate solutions are those computed by the third technique with which we deal, simulation. Because simulation of reservoir behaviour can be thought of as controlled sampling from the underlying stochastic process, it might be thought that more samples will, in aggregate, produce more accurate results. However, the statistics involved may be biased, even asymptotically, so it is useful to have some solutions (with known error bounds) with which to compare or test solutions obtained more conveniently by simulation or other approximate methods. We shall only touch on simulation: however, it must be stressed that where an input to a reservoir is described by the important class of mixed auto-regressive-moving-average (ARMA) models, the methods of difference or integral equations as used herein cannot yield solutions. (This is because the moving average c o m p o n e n t ensures that the model is not Markovian; indeed the current values of input are conditional on the entire history of the input in an ARMA model.) By contrast, simulation methods can handle ARMA inputs nearly as easily as independent inputs, with only a fractional increase in computational effort. Difference and integral methods can however handle the class of inputs described as auto-regressive, but with a huge increase in computational effort when compared to the independent input case. This cost increase demands that algorithms be found which can both reduce dimensionality and speed computation, so some attention is given to solving this problem.

271 PRELIMINARIES Before we proceed with a description of the algorithms and the results, we list some definitions and explanations for clarity of purpose.

Demand, draft. The demand 5 is assumed to be constant (all its variability being embodied in the input) and is the a m o u n t that one would like a reservoir to supply in a finite interval (t, t + 1]. The draft Yt, on the other hand, is a random variable which is conditional on the reservoir level and the input; Yt = 5 while there is storage available, but there will be occasions when the draft will be less than the demand. E l Y t] is thus less than 5 in a finite reservoir. N e t input. The net input Xt in an interval (t, t + 1] is defined as the gross input Zt less the demand 5. Thus E[Xt] = E[Zt] -- 5. (The term " n e t i n p u t " is borrowed from range and deficit analysis where E[Xt] = 5 invariably.) In a finite reservoir, X t is thus not truly a " n e t " input [Z t -- Y t ] , but since the adjective " a d j u s t e d " has come to acquire other connotations, and a thesaurus is not much help, we will use " n e t " . We use p = E [ X t ] , 02 = var[X t] and p = c o r r [ X t X t - 1 ] and assume stationarity throughout the sequel. We will derive results for three different marginal probability distributions for Xt; normal, log-normal and three-state discrete and treat both serially uncorrelated and correlated cases of each. Storage. Consider a reservoir of finite capacity V. Define: St+l = min[V, m a x ( 0 , S t + Xt)] Depending on the application, the initial storage So can be any point of the interval [0, V]. [We may be interested in the first passage time from 0 to V or from V to 0, or, more traditionally, in the equilibrium distribution of storage: lim P [ S t ~ s ] ,

(0~s~

V)

which is independent of So .]

Effective capacity. If, Z t and Yt o c c u r simultaneously, then the effective capacity of a reservoir of size V is equal to V. If, on the other hand, either Zt or Yt (or both) occur instantaneously, the former before the latter, then at t + 1 the level of the reservoir will be lower by ~ than it was the instant before the withdrawal commenced. The effective capacity of the reservoir is thus reduced by 5 in the case of "staggered" net input. All the results of the sequel will be interpreted in terms of an effective capacity equal to V. Thus they apply to a reservoir of actual capacity V when the net input is "simultaneous" (either instan-

272 taneous or evenly spread in (t, t + 1] ) and to a reservoir of actual capacity V + 5 when the net input is "staggered". The two results provide upper and lower bounds to what occurs in practice, where the net input is neither wholly staggered nor completely simultaneous.

Drift. It has been shown by Gomide (1975) that remarkable unification of the results of range and deficit analyses can be achieved by standardizing the variables. Mindful of this, we introduce a new term in reservoir studies, which we call the drift e, defined as the standardized mean net input p/a. It is thus the inverse of the coefficient of variation of the net input. Standard reservoir. We define a standard reservoir to be one of standardized capacity c = V/o fed by a standardized net input with drift e and unit variance. The scaling of the reservoir capacity b y a is in marked contrast to the c o m m o n practice of scaling the capacity by E[Zt] , the mean gross input. However, its introduction is justified because it effectively eliminates one of the independent variables from the computations with a concomitant reduction in cost and permits meaningful comparisons to be made of the behaviour of reservoirs fed by various input distributions. The probability of emptiness. In studies of reservoir reliability, a quantity of much interest is: v0 =

lira P [ S t = O ] t --+ o¢

the probability of emptiness. Also of interest is the quantity: Y0 =

lira P [ Y t ~ 5] t--~

the probability of failure (to supply the full demand). When comparing results from various types of input, it is important to know that the basis of comparison is sound. At early stages of this study it was thought that meaningful comparisons could be made by using Y0. However, this is not so, because in the case of the continuous input, v0 = Y0, whereas, by contrast, in the case of the discrete input, Y0 = Vo" P[Xt ~ 0] ~< v o (equality only holding when P[Xt ~ 0)] = 0). This suggests that when comparing the results from various input distributions one should use v 0 not y 0 as a basis. This is what is done in the sequel.

Mean first passage times. In studying the behaviour of a reservoir modelled by a finite Markov chain with discrete state--space and homogeneous (time invariant) transition probability matrix, we denote by mij the mean number of time steps taken to reach state i from state j for the first time, i,j = 0, 1 . . . . , c. In particular, if 0 is the e m p t y state of the reservoir/chain, moo

273

is the mean recurrence time of emptiness. Thus, if we have a discrete input to a suitably defined reservoir with a discrete state--space, then moo is none other than the reciprocal of v0, the probability of emptiness. Again, if c is the label of the full state, me0 is the mean first passage time from empty to full or the "mean time to fill", while conversely, m0c is the "mean time to e m p t y " , mcc is the mean recurrence time of being full, n o t usually of as much interest. The mean first passage time from full to e m p t y (or conversely) should be of some interest to hydrologists, because it is a statistic closely associated with the variable called the deficit by Gomide (1975), which in turn is the variable of interest in the "sequent peak algorithm" of mass-curve analysis. We are able to obtain approximate values of m0~ and m~0 for independent normal and log-normal inputs, and exact results for discrete inputs, but n o t for serially correlated inputs (except for the random walk) for reasons to be explained later where appropriate.

ALGORITHMS FOR DETERMINING THE DISTRIBUTION OF STORAGE

In this section are detailed the various algorithms for computing the equilibrium probability distribution of storage as a function of the standardized capacity c, together with the net input distribution parameters -- the drift e, the first serial correlation coefficient p and the skewness coefficient % Simulation

It has been shown elsewhere (Pegram, 1975) that the recurrence time of emptiness is approximately geometrically distributed in the independent input case. The quantity is n o t properly defined for a dependent input, but the following ideas carry over when p ¢ 0, nevertheless. Let too (the recurrence time of emptiness) be geometrically distributed with parameter p = moo = E [ t 0 0 ] . The waiting time to the kth emptiness has mean k p and variance k p 2 . The question one should ask before performing a reservoir simulation is: " H o w big would one choose k to ensure that the N trials sampled until the kth emptiness yield an estimate of p (fl = N / k ) , which has a given probability of lying within a given distance from p ? " When k is fairly large Feller (1968) shows that I] is normally distributed N ( p , lep: ) so t h a t the coefficient of variation of/5 is 1/X/~. If the coefficient of variation of ;fi is to be x%, then k = (100/x) 2 ; x -- 10% yields k = 100, x = 1% yields k -- 10,000. Of course, if one wished to be 95% sure that the value /5 is within x% of #, then k = (196/x) 2 , which almost quadruples the a m o u n t of computation. Now the expected number of trials to get three-figure accuracy by simulation in about two out of three cases (x = 1) is E [ N ] = E[klS] = p ( l O O / x ) : , so t h a t if p is a reasonable 100, E [ N ] = 106 , a respectable number for even

274

a fast computer. Thus, the prospect of ascertaining v0 (the probability of emptiness) as a function of c, e, p and 7 by simulation is a daunting one, to say the least, if one wishes much accuracy. We immediately find that subtle trends and variations due to changes in the independently chosen parameters are masked by the variability of/5. To counter this effect to some extent, a moderately long sequence of normally distributed pseudo-random numbers is generated and standardized. One then has a basic set with which to work, which is re-used for each simulation. (It must be stressed that this approach is only useful in comparative studies like this -- sensitivity of reservoirs to sampling variability should be handled differently, by taking white noise samples of the same length as the hydrologic sequence and feeding them through the simulation program to ascertain the effect of the original sample length on statistics of interest.) The simulation algorithm is of course trivial: given Xt one computes:

Ut+l = St + X t

(1)

If Ut+ 1 >~c, put St+ 1 = c ; i f

Ut+l <~0, put St+l = 0 ; otherwise St+l =

Ut+ 1 •

One records the number of times S is within any chosen sub-interval in (0, c) or at 0 or at c, to obtain a frequency distribution. The outcome of some simulations with n = 1000 and 7 = 0 appears in Table IX (p. 290).

Linear difference equations We start with the independent input case to introduce the ideas, then extend these to the dependent case. Using the storage equation (1), and putting g(u) as the p.d.f, of St, we can express the c.d.f. G(u) of St ÷ 1 as: c-

P[St+ 1 <~s] = P[St = 0 ] ' P [ X ' t ~ s ]

" P [ X t <~s --c]

+ P [ S t = c]

Again:

P[St+I < ~ s + h ]

+ fo ÷ g ( u ) ' P [ X t < ~ s - - u ] ' d u

= P[St = 0 ] ' P [ X t

for (0 ~< s ~< c)

(2)

<~s+h]

t " 12-

+ J0 ÷ g(u) • P [ X t <~s + h - - u ] • du

+ P[S t =c]'P[X t ~s +h--c] and taking differences, P[s ~St+

1

4 8 -~ hi

(--h % s < ~ c - - h )

= P [ S t -- 0] • P[s < X t <~s + hi °

+

÷g(u)'P[s--u
275 + P[St=c]

"P[s--c
]

(O~,:s~c--h)

If we divide the reservoir into k steps of equal size h, and put s -- ih, then for the interior "states": P[ih
~ ( i + 1)h]

1

= P [ S t = 0] " P [ i h < X t ~ ( i + 1)h] k-1

,(j+l)h

+ ~

(

g(u)'P[ih--u
1=0 ,2jh

+ P[St =c]" P[ih--c
<~(i + l ) t z - - c ] ,

i = 0, 1 , . . . ,

k --1

We approximate P [ j h ~ St + a <~ (J + 1)h] (which is an integral of a continuous distribution) by probability masses v i (t + 1) at the midpoints of the k interior sub-intervals of (0, c). Of course, P [ S t = 0] and P [ S t = c] are probability masses independent of the discretization, and their approximating equations follow directly from eq. 2. We can then write, following Moran (1959): v(t + 1) = Qv(t)

(3)

where v(t + 1) is the (k + 2) element vector probability mass function corresponding to the storage at time t + 1 and Q is a (k + 2)-square transition probability matrix whose elements are integrals of the input p.d.f. Thus if the c.d.f, of the input distribution is F(u), we find: F[O] qio

~---

i = 1, 2 , . . . ,

1 --F[c],

i= k + 1

for

=

i=1,2,...,k

1 --F[0],

i=k+l

l

F[(i--j+½)h]--F 1--F[c--(j+½)hl

k

i=0

F [ i h -- cl -- F [ ( i -- 1)h - - c l ,

F [ ( - - j + ½)h]

qij

i= 0

V [ i h ] - - F [ ( i -- 1)h],

F [-- c] qi,k+l

for

f!r (i--j--½)h],

i= 0 i=1,2,...,k i=k

+ 1

We find that reasonable accuracy can be obtained in most cases with k ~< 30; it is a matter of convenience as to which method is used to evaluate the equilibrium storage distribution vector: v =

lira v(t)

t ~o¢

276 One m e t h o d , the power m e t h o d , is to let v(O) be any vector whose (nonnegative) elements sum to uni t y and to c o m p u t e :

v(t + l) = Qv(t),

t = O, 1, 2, . . .

until there is no change. This usually takes a b o u t 20--30 iterations in the i n d e p e n d e n t input case, which is quite quick. Other m e t h o d s involve the solution o f te + 2 linear simultaneous equations, for the u n k n o w n v. We put Tv = e, where T is a matrix consisting of Q - I, with the last row replaced by a row o f ones, and e is a vector whose last element is one, the rest void. Gauss elimination or direct matrix inversion accomplish the solution conveniently. When the input is a lag-one auto-regressive time series, we have a multiple integration problem. The choice of the serially correlated log-normal distrib u t io n to describe the input as used in this paper was condi t i oned partly by the fact th at an algorithm for integrating the bivariate normal density was readily available. J ohns on and Kotz (1972) discuss m any algorithms, but the one chosen was tha t r e c o m m e n d e d by Sowden and Ashford (1969). If ¢(x, y, p) is the standardized bivariate normal density, then:

L(h, k, p) A fh

(p(x, y, p) • dy" dx { c2t --le

= f ] ~ Z ( t ) " c~P( ¢l ~, tl-' -ch~)

dt

where

Z(t) = (27r)-l/2exp( - t2/2) and t

• (t) = f

Z(x).dx

el = c2 -- X/P

ifp~0

and

el

= --¢2 = X/~---P if p < 0

The integration is p e r f o r m e d numerically by means of Gauss--Hermite quadrature, with 30 terms giving acceptable accuracy (6-figure) for I P I ~ 0.9. To define the algorithm for c om put i ng the marginal distribution (discretized) o f S t as t -~ o% we use the same discretization of the reservoir into k + 2 states as in the i n d e p e n d e n t case but instead of a univariate density function, work with the bivariate (discretized) density function of (St ÷l, St), which we call w(t + 1). (Note t hat in continuous reservoir studies, choosing the pair (St+ 1 ,St) rather than (St÷l ,Xt) as suggested by L l oyd (1963) for the discrete input, makes for a considerable reduction in the size of the matrix Q and the joint probability distribution vector.) The elements of w are:

wij = P [ ( i - - 1 ) h ~ S t + l

~ih,(j--1)h~St~jh],

i,j= 1, 2 , . . . , k

If i = 0 (or k + 1), the term "St+ 1 = 0 " (or c) appears in the bracket and

277

similarly for j and St. E v i d e n t l y we can r e c o v e r the marginal d i s t r i b u t i o n o f St b y s u m m i n g a p p r o p r i a t e l y : k+l

Vi =

L Wij j=O

T h e equilibrium d i s t r i b u t i o n o f (St+ 1 ,St), w m u s t satisfy Qw = w, w h e r e Q is a suitable t.p.m., the e l e m e n t s o f which are integrals o f the bivariate i n p u t d i s t r i b u t i o n ( n o r m a l or log-normal as the case m a y be). T h e integrations are d o n e over rectangular sub-regions o f the variable space (Xt+l , X t ) and this can be best e x p l a i n e d b y considering the m a s k (or t e m p l a t e ) f o r a small dimensional e x a m p l e . In Fig. 1 we have c h o s e n k = 4 and have labelled some o f the integrals tit c o r r e s p o n d i n g to the wit f o r a given starting state St-1 • D e p e n d i n g on w h e t h e r St-1 is 0, 1, 2 . . . . . 5, the m a s k is placed so t h a t the circles m a r k e d 0, 1, 2 , . . . , 5 c o i n c i d e with the (Xt+l,Xt) origin. R e t u r n i n g t o general k again, for each value o f St-1 (0, 1 , . . . , k + 1) the tit, i, j : 0, 1, . . . , k + 1, are evaluated by c o m p u t i n g L(h, k, p) using suitable values o f h, k and p. Once all the integrations have been p e r f o r m e d , f o r m t h e sums: k+l

t.j

=

Y i=0

T h e e l e m e n t s o f Q are: q i r I js

P[(i -- 1)h < St+ 1 ~ ih, (r -- 1)h < St <~rh I(J -- 1)h < St <~jh,

=

(s -- 1)h < S t_l <~sh]

t~5

t~

c

tzs

tso

St- i~-~,

tos 1

c

t3o

IAc: %

xt

4

tzo

to= too

to~

Fig. 1. T e m p l a t e f o r t h e e v a l u a t i o n o f t h e i n t e g r a l s tij o f t h e b i v a r i a t e d i s t r i b u t i o n o f ( X t , X t + 1 ) in t h e c a s e k = 4. Origin shown located at point 1 f o r t h e c a s e S t _ 1 ---- 1.

278

which are non-zero if r = j and zero otherwise, so in general we have qi~ (say) as non-zero elements. Now if we have chosen S t - i = s, then we compute qij~ = tij/t • j, so t ha t each column of Q is normalized, and Q is a bona fide t.p.m. T o summarize (working in discrete variables, and avoiding the boundaries for simplicity) : w O = P[St+ 1 = i , S t = j ]

8

=

P[St+ 1 = i , S t = j , S t _ 1 = s] = ~ P[St+ 1 = i [S t = j , S t - 1 = s]" P [ S t = j , S t - 1 8 P[X t=i-j[Xt_

1 =j--s]'P[S

8

t=j,St-

l=s]

=s]=

= ~. q i i s ' W i s 8

as St+ 1 = St + Xt e x c e p t at the boundaries. The qijs (elements of Q) are thus conditional probabilities ( X t [ X t - 1 ) which we have interpreted as normalized integrals of the bivariate continuous distribution, ensuring that the partition of the ( X t , X t - x ) plane is done in the mo s t economically accurate way possible. As to solving the system for w, recall that k ~ 30 for the univariate case. T o get comparable accuracy, the vector w would have to consist of ~ 1 0 2 4 elements, so th a t Q would contain a bout 33,000 non-zero elements. The p o w e r m e t h o d was used to solve for w. Only the positive elements qijs of Q were stored together with w(t + 1) and w(t), and then w(t + 1) was comp u t e d from Qw(t) with as much e c o n o m y as could be made. To test for convergence, the ratio:

max [wii(t + 1)/wij(t)]/min [wii(t + 1)/wiy(t)] ij

ij

was c o m p u t e d after every 10th multiplication. When it was less than 1 + t in magnitude (with t = 10 -7 ), convergence was assumed. Depending on the degree o f departure from normality and the magnitude of c and p, the number o f iterations varied from 20 to about 70. If a solution t o o k m ore than 100 iterations to converge, it was f o u n d t hat the particular problem was suffering from instability caused by t oo coarse a discretization. The result was th at only certain regions of the parameter space were able to be covered w i t h o u t incurring very large c o m p u t e r costs. The results of the c o m p u t a t i o n s are listed after the n e x t section.

Integral equations Starting with the i n d e p e n d e n t input case, we rewrite eq. 2, putting f(x) and F(x) as the p.d.f, and c.d.f, of the input Xt (not necessarily normal) and g(s) and G(s) as the p.d.f, and c.d.f, of the storage St :

279 ff c

G(s) = G(O) • F ( s ) + J,÷ g ( u ) F ( s -- u) • d u + [1 -- G(c-)] F ( s -- c),

(0 ~ s ~ c)

(4)

Integrating the r.h.s, by parts we get: C-

a(S)

= F ( s -- c) + Jo f(s -- u) " a ( u ) " du

(5)

which is a Fredholm integral equation of the second kind. Some exact solutions to eq. 5 exist in special cases, but a more flexible approach is the numerical solution of the integral equation. The literature suggests that solutions be obtained by quadrature, by expansion methods or by successive approximations. All three methods were tried, and quadrature was found to be the neatest, fastest and most easily implemented in this context, the reason being the lack of an explicit expression for the input c.d.f. (An exact storage distribution was obtained for an exponential input, but only because the exponential function is so easy to handle.} Dividing the interval (0, c) into k equal slices, we evaluate an approximation H ( s ) to G(s) at each of the boundaries of the sub-intervals. Thus if h = c / k , we evaluate the approximation at s = ih, i = 0, 1, . . . , k. Using a numerical quadrature formula, such as trapezoidal or Simpson's, the integral on the right of the eq. 5 is evaluated in terms of the unknowns H ( ] h ) . Thus we get k + 1 linear simultaneous equations in the k + 1 unknowns. If the weights of the quadrature formula are wj, j = 0, 1 , . • . , k, then we write: k

H ( i h ) = F ( i h -- c) + ~

wif[(i --j)h] "H(jh)

.i=O or

x = b+KWx

where x i = H(ih);

bi = F(ih--c);

Ki~ = f [ ( i - - j ) h ]

and )

W = diag(wi) Then x

(6)

= [I--KW]-lb

is the solution for the approximating c.d.f. H [ s ] . The advantage of the quadrature m e t h o d is that one immediately has an interpolation formula which yields H ( i h ) at the nodes, and employs all the information in the kernel to interpolate between them, i.e.: k

H(s) = F(s--c)+

~. w j ' f [ s - - j h ] ' H ( j h ) i=0

280 [This r e m a r k a b l y e f f i c i e n t m e t h o d can easily be m o d i f i e d t o include a variable d r a f t rule, e v a p o r a t i o n , etc. (Kartvelishvili, 1969, p. 115). It is relatively s t r a i g h t f o r w a r d to e x t e n d to the d e p e n d e n t n o r m a l , b u t the d e p e n d e n t logn o r m a l i n p u t appears t o require a lot m o r e e f f o r t . ] T o gain a c c u r a c y , the d e f e r r e d a p p r o a c h t o the limit ( h 4 - e x t r a p o l a t i o n ) was e m p l o y e d . This necessitates the successive d o u b l i n g of k (halving the interval), and yielding a sequence o f estimates of v0, the p r o b a b i l i t y o f failure: T(1) . . . Tim) which usually converges f r o m either above or below. Baker ( 1 9 7 7 ) suggests the c o m p u t a t i o n of:

U(i) = 2T(i) -- T(i - - 1 ) ,

i=2,3

....

,m

as giving a sequence which converges f r o m the o t h e r side. These are t h e n u p p e r and lower b o u n d s t o the desired value. In o u r c o m p u t a t i o n s , m = 4, and we used k - 2, 4, 8, 16 in t u r n t o o b t a i n our T(i). T h e d e f e r r e d app r o a c h t o the limit in c o m b i n a t i o n with S i m p s o n ' s yields: T¢°)(i)

= T(i),

T(1)(i)

=

i-- 1,2,3,4

[16T¢°~(i) -- TC°~(i --1)] /15,

T¢2~(i) = [ 6 4 T ( l ~ ( i ) - - T ¢ I ~ ( i - - 1 ) ] / 6 3 , T¢3~(4) =

i=2,3,4 i=3,4

[256T¢2~(4)--Tc2~(3)]/255

We find t h a t T~3)(4), our closest estimate of v0, lies b e t w e e n U(4) and T(4), which give us b o u n d s on the a c c u r a c y o f our solution. T h e results for a few c o m p u t a t i o n s are r e p o r t e d in the n e x t section b u t one, where we c o m p a r e t h e m with results c o m p u t e d by simulation and diff e r e n c e e q u a t i o n s ; h o w e v e r , an e x a m p l e o f a h a n d calculation m a y help the ideas. F o r a n o r m a l i n p u t with e = 1, p = 0 fed into a s t a n d a r d reservoir with c = 1, we get, f o r k = 2: 0.2420

0.3521

0.3989

K = 0.1295

0.2420

0.3521

0.0540

0.1295

0.2420

)

W = diag[1/6

4/6

b

0.0668

0.0228] w

0.0868

0.1889] T

= [0.1587

1/6]

yielding x

= [0.0333

f r o m eq. 6. T h u s the p r o b a b i l i t y o f emptiness is e s t i m a t e d t o be 0 . 0 3 3 3 in this case ( c o m p u t e d b y hand-held calculator), which c o m p a r e s well with the " a c c u r a t e " result 0 . 0 3 3 2 5 6 ! B e f o r e we go o n to presenting the results, we give a q u i c k outline o f the

281 a l g o r i t h m s f o r o b t a i n i n g m e a n first passage times. T h e m a t e r i a l exists elsew h e r e , b u t it is c o n v e n i e n t to have all the ideas " u n d e r o n e r o o f " as it were.

ALGORITHM FOR EVALUATION OF MEAN FIRST PASSAGE TIMES OF A MARKOV CHAIN K e m e n y and Snell ( 1 9 6 0 ) detail the d e v e l o p m e n t of the algorithm. We q u o t e t h e i r results a n d s h o w h o w t h e y can be applied in storage analysis. L e t a discrete-state s t o c h a s t i c processs { S t } be a lag-one finite M a r k o v chain w i t h n + 1 states, d e f i n e d b y its (k + 2)-square t r a n s i t i o n p r o b a b i l i t y m a t r i x Q. {S,~} will be a s t a t i o n a r y process if Q is c o n s t a n t . T h e m a r g i n a l p r o b a b i l i t y d i s t r i b u t i o n v e c t o r of S t , d e n o t e d b y p ( t ) = { P [ S t = i] } i = 0, 1, 2, . . . , k + 1, is given b y p ( t ) = Q p ( t - 1) = Q t p ( O ) w h e r e p ( 0 ) is the initial distribution vector. I f Q has at least o n e positive diagonal e l e m e n t a n d n o zero r o w s or c o l u m n s , t h e n S t is ergodic, w h i c h implies t h a t p ( t ) - ~ 7r, (the e q u i l i b r i u m d i s t r i b u t i o n v e c t o r o f S) as t --> ~ . I t follows t h a t QTr = 7r, w h i c h w h e n c o m b i n e d with the c o n d i t i o n l'Tr = 1 ( w h e r e 1' is a r o w o f ones), yields ~ on s o l u t i o n o f the set o f k + 2 i n d e p e n d e n t linear s i m u l t a n e o u s e q u a t i o n s , as we saw earlier. K e m e n y a n d Snell ( 1 9 6 0 ) give the f o l l o w i n g a l g o r i t h m f o r e v a l u a t i o n o f t h e m e a n first passage t i m e m a t r i x M = { m i j }, i , j = O, 1, . . . , k + 1, w h e r e each e l e m e n t mij is e q u a l to the m e a n n u m b e r of steps t a k e n t o go f r o m j t o i f o r t h e first time. ( I t is usual t o call m i i t h e m e a n first r e c u r r e n c e t i m e f o r s t a t e i, b u t we will dispense with t h a t d i s t i n c t i o n w h e r e we are talking collectively o f m i j . ) T h e y find M as follows. D e f i n e the " f u n d a m e n t a l m a t r i x " Z as: Z =

[I-Q+~rl']-I

then M

= D[I--Z

(7)

+ dl']

w h e r e D is a diagonal m a t r i x w i t h 1/Tr i as the ith e n t r y a n d d is a v e c t o r f o r m e d f r o m the d i a g o n a l e l e m e n t s o f Z. When { S t } is a r a n d o m w a l k w i t h p a r t i a l l y r e f l e c t i n g barriers at 0 a n d k + 1, so t h a t : P[St =liSt-1 P[St=O[St-

=0] ! =0]

P [ S t = i + 1 [ S t = i]

= P[St=k+ -- P [ S t = k [ S t - 1 = p

and

l[St-1

=k+

=k+l]

l]

= p

= q -- 1 - - p

P [ S t = i -- 1[ S t_l = i]

= q

f o r i -- 1, 2 , . . . ,

k

t h e n { S t } is a finite discrete ergodic M a r k o v chain. K e m e n y a n d Snell ( 1 9 6 0 ) give the f o l l o w i n g e x p l i c i t f o r m u l a e f o r mij f o r this case.

282 I f p = 1/2:

i=j

k+2, rni~ = l ( 2 k - - i + 3 ) i - - ( 2 k - - j + 3 ) J '

i>j

[j(j + 1) -- i(i + 1),

(8)

i
I f p :¢: 1/2: r k+2 -- 1 1

r~ =I

mi,

11

[ P--q[

r-) ,

i= j

[r k - j + 2 - r k - ' + 2 r---1 . . . .

j--i)

]

(i--j)],

(r_ 1)r~+j],

i>j

(9)

i
where r -- p/q. T h e s e f o r m u l a e will p r o v e useful w h e n we s t u d y the discrete finite reserv0ir fed b y a + 1 input. In the case o f the i n d e p e n d e n t continuous net i n p u t ( n o r m a l or logn o r m a l ) we have already indicated in the previous section h o w the approxim a t e t.p.m. Q is c o n s t r u c t e d b y discretization. F o r m a l l y , the c o n t i n u o u s and discrete cases are t r e a t e d the same; w h e n c o m p u t i n g M. We e m p l o y eq. 7 t o c o m p u t e moo, m0c, me0 and mc~ numerically, increasing k until we are c o m f o r t a b l e a b o u t the accuracy. Here the subscript c d e n o t e s the full state in a reservoir fed b y a c o n t i n u o u s input. When we consider a serially d e p e n d e n t i n p u t , we c a n n o t define m e a n first passage times f r o m one group o f states to a n o t h e r unless these are " l u m p a b l e " . [See K e m e n y and Snell ( 1 9 6 0 ) f o r d e f i n i t i o n o f and restrictions on l u m p a b i l i t y . ] In o n l y one case which we t o u c h in this p a p e r d o we have l u m p a b l e states at the boundaries, and t h a t is in the case o f the c o r r e l a t e d r a n d o m walk (+ 1 input). This is because the c o l u m n s o f Q c o r r e s p o n d i n g t o the states (St+l,St)= ( 0 , 0 ) and ( 0 , 1 ) are the same. So are the c o l u m n s (St+l,St) = (k + 1 , k ) and (k + 1 , k + 1). We can thus d e f i n e the states St+l = 0 and St+l = k + I as the ( l u m p e d ) b o u n d a r y states o f the bivariate chain (St+l,St). It is o n l y in this case t h a t we can so do ( e x c e p t for the d e g e n e r a t e case p = 0 w h e n all pairs o f c o l u m n s are the same so t h a t the chain reduces t o a univariate chain). We shall give a fuller t r e a t m e n t o f the discrete i n p u t below. RESULTS T o o b t a i n a q u a n t i t a t i v e idea o f the a c c u r a c y o f the a p p r o a c h e s , we first r e p o r t the result o f c o m p u t i n g the p r o b a b i l i t y o f failure v 0 o f a s t a n d a r d

283 reservoir o f size c fed b y an i n d e p e n d e n t n o r m a l i n p u t with drift e. T h e comp u t a t i o n s are d o n e b y n u m e r i c a l l y solving the relevant d i f f e r e n c e and integral equations. T h e reservoir is successively divided into k = 2, 4, 8, 16 slices, and we c o m p u t e u p p e r and l o w e r b o u n d s t o v 0 as a f u n c t i o n o f k. O n l y t h o s e sequences which are m o n o t o n i c are r e p o r t e d , because the discretiza t i o n and r o u n d o f f in some cases gives results which do n o t yield m o n o t o n i c sequences and these are n o good for calculating b o u n d s . Taking the m o r e a c c u r a t e o f the t w o sets o f results, an estimate o f v o is given in each case. T h e results a p p e a r in Table I. One o f the o u t c o m e s o f this c o m p a r i s o n is t h a t we can e x p e c t " a c c u r a t e " results for small c and " a p p r o x i m a t e " results for c ~ 2, w h e n using the finite-difference algorithm. T h e remaining results will all be r e p o r t e d as m e a n first passage times (to three significant figures). Where we have been u n a b l e to attain this a c c u r a c y , the result will be o m i t t e d , and w h e r e it is o n l y within 10%, we will insert the figure, b u t p u t it in parentheses. T h e r e are t h r e e main groups, c o m p u t e d b y d i f f e r e n c e equations. Firstly, we have the results for i n d e p e n d e n t n o r m a l n e t input, giving the m e a n first passage times moo, moc, mco and m~c as f u n c t i o n s of the standard c a p a c i t y c and drift e (Tables II--V). T h e second group will give moo as a f u n c t i o n o f c, e and ~ for an i n d e p e n d e n t log-normal input. F o r this case, we n e e d t o find the three p a r a m e t e r s o f the log-normal distribution, n a m e l y fl, p and 02 , in terms o f E[Xt] = O, v a r [ X t ] = 1 and % h e n c e we need the positive r o o t o f the cubic e q u a t i o n : 72(z-1)

3 = (z 3 - 3 z + 2 )

2

Having solved f o r z, we get: 02 = l n z ;

p =

--[ln(z z-z)]/2;

fl = - - e x p ( o 2 / 2 + p )

The d i s t r i b u t i o n of X t is shifted by e to a c c o m m o d a t e n o n - z e r o d r i f t (Table VI). T h e third group o f results r e p o r t s moo as a f u n c t i o n o f c, e, 7 and p for a serially c o r r e l a t e d log-normal input, in which the c o r r e l a t e d n o r m a l appears as a special case with 7 = 0 (Tables VII--IX). Finally, Table X gives a comparison b e t w e e n simulation and d i f f e r e n c e e q u a t i o n results for interest's sake.

Independent, normally distributed net inputs (Tables II--V} N o t e t h a t , because o f the s y m m e t r y o f the n o r m a l d i s t r i b u t i o n : [m00[e]

= [mccl--e]

and

[mco[el

= [mocl--e]

Independent, log-normally distributed net inputs T h e results are s h o w n in Table VI.

284

"0 e~

~.

L~

~

O0

O0

~

~"

~

~0

"~

~ 0

0 0

Cq 0

0 0

0 0

"0

0 0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

~ 0 0 0

0 0 0 ~

0 0 0 0

0 0 ~ 0

3 "0

0 0 0 0

3 0

e~ Cq

e~ 0

0 Z~

0

e~ 0

O

LO O

0 ~--~

C)

L~

~

r.,)

._

0

C~q

LO 0

O. ,-t

285

0

b0 0

~

~

QO b - ~0

0 0 0

~.qq

0 ~D bCq 0

o

~,

0 0

0 cO

,

,,

0 0 0

6

b-

<:~

0 0 0

0 0 0

b00

0 0 0

q 0

00

0

T-~

286 T A B L E II moo, mean recurrence time of emptiness

c• 0.25 0.5 1 2 4 8 16

0

0.2

0.4

0.6

0.8

1.0

2.22 2.48 3.08 4.47 7.31 13.0 24.5

2.74 3.19 4.38 8.00 21.6 115 2,630

3.49 4.27 6.60 16.1 84.7 1,920 95 • 104

4.60 5.95 10.5 35.7 397 41,700 45 • 107

6.38 8.61 17.4 84.6 2,060 101 • 104 --

8.85 12.9 30.1 210 11,300 26 • 106 --

T A B L E III

moc, m e a n first p a s s a g e t i m e f r o m full to e m p t y

0.25 0.5 1 2 4 8 16

2.46 3.07 4.75 10.0 26.6 83.6 288

0.2

0.4

0.6

0.8

1.0

3.02 3.89 6.47 15.9 59.9 408 10,300

3.83 5.11 9.22 28.0 176 4,250 21 " 1 0 s

4.99 6.94 13.8 54.2 653 69,400 76 "107

6.72 9.76 21.6 114 2,870 14 ° 10 "~ --

9.35 14.2 35.3 258 14,100 33 • 106 --

T A B L E IV

moo, m e a n first p a s s a g e t i m e f r o m e m p t y to full

c• 0.25 0.5 1 2 4 8 16

0

0.2

0.4

0.6

0.8

1.0

2.46 3.07 4.75 10.0 26.6 83.6 288

2.06 2.49 3.64 6.86 14.9 33.5 73.2

1.77 2.09 2.91 5.06 9.87 19.8 39.8

1.56 1.79 2.40 3.96 7.27 13.9 27.3

1.40 1.58 2.04 3.24 5.74 10.7 20.7

1.29 1.42 1.78 2.74 4.75 8.75 16.7

287 TABLE V mcc, m e a n first r e c u r r e n c e t i m e o f fullness

c•

0.25 0.5 1 2 4 8 16

0

0.2

0.4

0.6

0.8

1.0

2.22 2.48 3.08 4.47 7.31 13.0 24.5

1.86 2.00 2.31 2.85 3.50 3.94 4.07

1.60 1.68 1.84 2.05 2.20 2.26 2.27

1.42 1.47 1.54 1.63 1.67 1.68 1.68

1.29 1.32 1.36 1.39 1.40 1.41 1.41

1.20 1.21 1.23 1.24 1.25 1.25 1.25

Serially correlated log-normal net input (Tables VII and VIII) The algorithm does n o t behave well when 7 is large. The figures in Table VII were all c o m p u t e d by the power m e t h o d applied to the difference equations. It will be noticed t hat some of the figures for p = 0 do n o t agree in the last digit with values in Table VI. This is because those in Table VII were c o m p u t e d from a bivariate input with p = 0, com bi ned with the fact t h a t different m e t h o d s were used for each calculation. The results for c larger than 2 are inaccurate (to 3 significant figures) because it was n o t economical to make k large enough. Hence they will n o t be report ed here.

Results o f simulation, normal net input (Table IX) A standard reservoir of capacity c = 2 was fed with several sequences of auto-regressive samples of length 1000, all generated from the same set of pseudo-random normal variables. These values appear in Table IX, together with the relevant figures from Table VIII in brackets for comparison. Note t h a t the simulation figures are not accurate, but reflect the true behaviour quite well for small e; the consistent over-estimation of moo could be ascribed to some sample correlation in the standardized white noise sequence, in addition to sample bias.

EXACT RESULTS RESERVOIR

FOR

DISCRETE

A P P R O X I M A T I O N S TO THE

CONTINUOUS

To apply the t h e o r y of r a ndom walks to storage analysis, we must take care with the specification of the state--space of the discrete reservoir. To have the correct properties, the discrete reservoir of "size" k must have k + 2 states, 0, 1, . . . , k + 1 where 0 and k + 1 are the " e m p t y " and " f u l l " states. (This is justified if it is realized that for a symmetric i nput t o a continuous reservoir, v0 -~ 0.5 as c ~ 0. The equivalent discrete reservoir is one

288

% ! I

Itll

0o0 O 0

I I

0

~o

~

%

%

%

%

~D

%

%

%

00 0 ~

O,1000

.r-.~

~4D 01 00O,1

~U %

%

%

00 Q

0'~

O ~

C

°~

%

%

% O ~ O ~

~0 ~.. O0 0

0.. o

% o ~J

0

~

0

~ 0 . o .

0

~ o

,~,-.~ 00 ~

o~o

~

0

~ 0 ~

0 0

0

0

0

289

T A B L E VII moo, m e a n r e c u r r e n c e t i m e o f e m p t i n e s s as a f u n c t i o n o f e, 3' and p for a s t a n d a r d reservoir o f size c = 1 e

7

P 0

0.2

0.4

0.6

0.8

0

0 0.5 1 2 4

3.08 2.83 2.65 2.45 2.29

2.98 2.75 2.58 2.39 2.22

2.89 2.68 2.52 2.33 2.15

2.80 2.61 2.46 2.26 2.07

2.70 2.51 2.37 2.17 1.99

0.2

0 0.5 1 2 4

4.38 4.11 3.98 4.04 4.40

4.09 3.84 3.70 3.65 3.83

3.86 3.62 3.47 3.36 3.41

3.64 3.42 3.26 3.10 3.07

3.42 3.20 3.04 2.85 2.76

0.4

0 0.5 1 2

6.59 6.64 7.13 9.62

5.93 5.87 6.09 7.34

5.40 5.27 5.33 5.94

4.94 4.77 4.73 4.98

4.49 4.28 4.18 4.21

0.6

0 0.5 1 2

10.5 12.1 16.2 44.9

9.06 9.94 12.1 24.3

7.93 8.37 9.53 15.3

6.98 7.15 7.72 10.5

6.12 6.08 6.29 7.59

0.8

0 0.5 1 2

17.4 24.9 50.6 (800)

14.6 19.0 31.2 (200)

12.2 14.8 20.9 (80)

10.3 11.7 14.8 (40)

8.67 9.27 10.7 (20)

1.0

0 0.5 1 2

30.1 58.7 239 (200"103 )

24.5 41.2 113 (15"103 )

19.8 29.3 60.0 (2,000)

16.0 21.2 35.0 (400)

12.8 15.3 21.4 (100)

with k = 0, so t ha t k ÷ 2 = 2 and v0 -- ½.) Transitions from one state to a n o t h e r will only occur if the i nput consists of am ount s which are multiples o f the difference between adjacent storage levels. Thus, the " d i s t a n c e " between the adjacent states of a discrete reservoir are all equal. By contrast, in the finite-difference a ppr oxi m a t i on to the c o n t i n u o u s reservoir, the distance between the adjacent interior states is Ac, while the distance from b o u n d a r y to adjacent interior states is Ac/2. Despite this disparity, the results for the c o n tin u o u s and discrete reservoirs compare favourably. In the following, we must distinguish between cases where e -- 0 and e =~ 0.

290 T A B L E VIII moo for c = 2 e

~

p 0

0.2

0.4

0.6

0.8

0

0 0.5 1 2

4.48 4.13 3.87 3.52

4.07 3.76 3.52 3.21

3.74 3.45 3.23 2.96

3.46 3.19 2.99 2.73

3.17 2.93 2.75 2.50

0.2

0 0.5 1 2

7.99 7.81 7.79 8.20

6.55 6.27 6.14 6.26

5.57 5.26 5.08 5.03

4.82 4.53 4.33 4.17

4.19 3.91 3.71 3.49

0.4

0 0.5 1 2

16.1 18.4 22.5 40.3

11.6 12.3 13.8 20.5

8.90 9.01 9.56 12.2

7.12 6.99 7.11 8.12

5.76 5.53 5.45 5.68

0.6

0 0.5 1 2

35.6 55.3 104 (710)

22.3 28.9 43.6 (170)

15.3 17.7 23.0 (60)

11.1 11.9 13.9 (25)

8.28 8.40 9.00 (12)

0.8

0 0.5 1

84.1 211 864

46.2 82.8 215

28.0 40.6 76.5

18.3 22.9 34.2

12.4 14.0 17.4

1.0

0 0.5 1

54.8 111 395

32.1 50.6 llG

19.6 25.6 41.7

208 1,030 15,200

102 293 1,890

T A B L E IX C o m p a r i s o n o f s i m u l a t i o n a n d d i f f e r e n c e e q u a t i o n results (in b r a c k e t s ) for c -- 2 e

p 0

0.2

0.4

0.6

0.8

0

4.63 (4.48)

4.24 (4.07)

3.97 (3.74)

3.47 (3.46)

3.08 (3.17)

0.2

9.26 (7.99)

6.99 (6.55)

5.99 (5.57)

5.13 (4.82)

4.27 (4.19)

0.6 1.0

66.7 (35.6) 500 (208)

40.0 (22.3) 333 (102)

21.3 (15.3)

14.5 (11.1)

200 (55)

83 (32)

10.2 (8.28) 32.3 (20)

291

I n d e p e n d e n t discrete + 1 i n p u t When discussing first passage times, we will c o n t i n u e t o use t h e s u b s c r i p t c t o d e n o t e full; t h e n f r o m eq. 8 w i t h k = c, a n d e = 0: moo = mcc = k + 2

and

m0c = me0 = ( k + l ) ( k + 2 )

(10),(11)

N o t e t h a t a l t h o u g h a discrete reservoir c a n n o t be d e f i n e d f o r n o n - i n t e g e r k, if we allow t h e a b o v e e x p r e s s i o n s to be i n t e r p r e t e d f o r 0 < k < 1, t h e y closely a p p r o x i m a t e the results f o r t h e i n d e p e n d e n t n o r m a l d i s t r i b u t i o n input. T h e p r o b a b i l i t y d i s t r i b u t i o n o f t h e discrete -+1 i n p u t w i t h n o n - z e r o d r i f t is given b y : P[Xt=+I]

= p,

P[Xt=--I]

= q = 1--p

withp:/:l/2

T o c o m p a r e the results of the n o r m a l a n d discrete inputs, we m u s t define the s t a n d a r d r e s e r v o i r f o r discrete i n p u t s w h e n e ve 0. Firstly, t~ = P - - q, a = 2x/'P-q, so t h e drift: e = (p -- q)/(2x/rPq -) S e c o n d l y , f o r a given k, the s t a n d a r d c a p a c i t y c = k / o a n d vice versa. When e ¢ 0, we have r = p / q =/=1. In t e r m s o f e: p = ( v z i + e 2 + e ) / ( 2 x / ~ - + e 2 ); (p -

q)

=

e/x/i

r = (x/l- + e 2 + e)/(X/~ + e 2

-

-

e)

+ e2

so t h a t we get, f r o m eq. 9: moo

=

(r/+2 --1)/(r--I)

moc

=

[m00--(k

rncc

= moo~ r k + l

mco

= [k + 2 - - m c ~ ] / ( p - - q )

+2)]/(p--q) p g=q

T h e s e e x p r e s s i o n s clearly i n t e r p o l a t e f o r n o n - i n t e g e r k so t h a t we can use t h e m to c o m p u t e first passage t i m e s f o r given values o f c a n d e, b y w h i c h t h e y are c o m p l e t e l y specified. I t m a y be m o r e c o n v e n i e n t t o c o m p u t e c f r o m a given moo a n d e. In t h a t case:

a = 1/X/~+e

2

and

k = ln[l+m00(r--1)]/lnr--2,

w h e n c e c = k / a is t h e s t a n d a r d c a p a c i t y . Discrete three-state i n p u t s o t h e r than + 1 O f course, we are n o t r e s t r i c t e d t o the -+1 process. I f we c h o o s e X t t o h a v e a b i n o m i a l (2, p ) d i s t r i b u t i o n , t h e n the f o l l o w i n g r e l a t i o n s h i p s hold.

292 o = "v/-2pq

p = (1 + e / X / ~ + e2)12

e = (p--q)/o

q = (1--e/V~+e2)/2

r

= [ ( V ~ + e 2 + e ) / ( x / ~ + e 2 -- e)] 2

= p2/q2

For example, when c = 4 and e = 0.2, we get a = 0.70014, k = 40 = 2.80056, r = 1.75737, yielding moo = 18.46. Comparing this with 19.51 for the +1 input, and 21.60 for the normal, there is n o t a large disparity. By contrast, for c = 2 and e = 1, the normal yields moo = 209.2, the binomial 314.1 and the +1 input gives moo = 84.9, so the expressions appear n o t to be so close when e (and hence moo ) is large. Serially c o r r e l a t e d d i s c r e t e i n p u t s w i t h arbitrary d r i f t

Since the a p p r o x i m a t i o n of the continuous reservoir by the discrete model is fairly satisfactory in the case of i n d e p e n d e n t inputs, some confidence is felt th at the same will occur when the inputs are no longer independent. For the discrete input, the analogue of the lag-one auto-regressive model is the lag-one Markov chain, specified by the transition probability matrix: L = p I + (1 - - p ) 7rl'

(12)

where p is the first serial correlation coefficient; I is the identity matrix; and is the equilibrium ve c t or probability distribution of { X t }. This input model was i n t r o d u c e d by Pegram (1974) where it was applied to reservoir problems. Explicit expressions for Y0, the probability of failure to m eet full demand, were given by Pegram (1978), however, as we have seen above, we need to c o m p u t e v0 (or its inverse m00) rather than Y0A brief outline of the derivation of the expressions will be given, and we then explore the behaviour of .moo as a f u nct i on of c, e and p. The joint equilibrium vector probability distribution of [ S t , X t ] is v = {vii } where: vi~- =

lim P [ S t = i, X t = j ],

i = O, 1, 2 , . . . ,

k + 1;

j = -- 1, O, + 1

t-+co

(Note that when X t is a three-state input it is m ore economical to work with the joint distribution of ( S t , X t ) rather than (St+l , S t ) when k ~ 1.) Then the equilibrium storage distribution is given by the elements: +1

vi =

lim

P[S t =

i]

=

t -+ ~

~

vii

j= -1

It was shown by Pegram (1974) t hat in particular: vo =

a/[a+

b(1 + h + h 2

+...+Xk-l)+dXk]

(13)

for a (k + 2)-state reservoir, where: a = 1/[p(1 --p)] b = (l+p)/(q+pp)

d = 1 / [ q ( 1 - - p)] ~ = (p +pq)/(pp

(14) +q)

293

Pegram ( 1 9 7 8 ) derived these expressions (14) w h e n the i n p u t has the particular f o r m given b y eq. 12. In t h a t case, ~ = [q (1 - - p - q ) p ] T is the equilibrium vector distribution of Xt. T h e r e are t w o cases o f interest, zero and n o n - z e r o drift. F o r zero drift, p = q and the s u b s t i t u t i o n o f eq. 14 i n t o eq. 13 yields, a f t e r simplification: moo

= mcc = 1 / %

(15)

= k + 2--pk

which in turn yields eq. 1O w h e n p = 0. N o t e t h a t this expression is quite i n d e p e n d e n t o f p or q and h e n c e the variance of the i n p u t . When e :~ 0, we find t h a t : moo

= (pX~+l/q--1)(pp

(16)

+ q)/(p --q)

which b e c o m e s (r k+2 - - 1 ) / ( r - - 1 ) , as e x p e c t e d w h e n p = 0, since in t h a t case, ?t = p / q = r. We also have t h a t (see eq. 9): mcc = m o o / ( p X k / q ) b e c o m e s m o o / r k+l w h e n p = 0

We n o t e t h a t for a three-state i n p u t p and q can in general be c h o s e n indep e n d e n t l y (subject to the c o n s t r a i n t t h a t all e l e m e n t s o f 7r are non-negative). H o w e v e r , if we specifically fix p + q = 1, t h e n we o n l y have one p a r a m e t e r , p, describing the equilibrium marginal d i s t r i b u t i o n o f the discrete -+1 input, while p fixes the serial correlation. We t h e n have t h a t , as b e f o r e : p = (x/l+e

2+e)/(2V~+e

2)

and

k = 2cV~

so t h a t moo is f o u n d in t e r m s o f c, e and p only. J u s t f o r fun, we indulged in some e x p l o r a t o r y c o m p u t a t i o n which y i e l d e d an expression for moc w h e n p =~ 0, e = 0, and the i n p u t is discrete -+1. It appears in Table X. In t h a t case the states St = 0 and St = k + 1 can be l u m p e d , so t h a t the m e a n first passage times are p r o p e r l y defined. T o find the f o r m u l a allgebraieally w o u l d have been tedious; it has been t e s t e d b y

TABLE X Mean first passage times o f a discrete reservoir fed by a c o r r e l a t e d Bernoulli i n p u t

p=O

p#=O

For e = O : moo

rnoc

k + 2 (k + 1)(k + 2)

k+2--pk [k/(1 + p ) + 11(1 - - p)] [k + 2 --pie]

For c:¢: 0: moo

(r h+2 - - 1 ) / ( r - - 1)

moc

[moo--(k+

2)]/(p--q)

(rX k+l -- 1)(pr -I- 1)l(r - - 1) use eq. 7

294

c o m p u t a t i o n for various values o f k, e and p, and was f o u n d t o be e x a c t in all cases tested. It also yields the c o r r e c t f o r m u l a w h e n p is p u t equal t o zero. T o s u m m a r i z e then, w e c o l l e c t the entries in Table X w h e r e k, p, q and r are d e f i n e d in t e r m s o f c, e, p f o l l o w i n g eq. 11 and k is d e f i n e d in eq. 14. Values f o r e =/: 0, p =/= 0 can be evaluated using eq. 7 and an appropriate Q, b u t since w e have n o figures with w h i c h to c o m p a r e t h e results o f such a c o m p u t a t i o n , it will n o t be u n d e r t a k e n here. S u f f i c e it t o say that in d o i n g s p o t c a l c u l a t i o n s c o m p a r i n g the f o r m u l a e o f Table X w i t h results for t h e n o r m a l input, w e find that t h e general c o m p a r i s o n is favourable, and in m o s t cases t h e discrete reservoir yields quite d e c e n t a p p r o x i m a t i o n s .

0

3 C

Fig. 2. The relationship between size (c) and drift (c) of a standard reservoir fed by independent discrete -+1 input (full lines) and independent normal input (dashed lines) as a

function of the mean recurrence time of emptiness (moo).

0

2

4

,6

32

6,

,26 2~6 ,,2 :

.

.

.

,o2, ~o,8 moo .

p.-o.

0.2

04 E 06

V,

\

\

"zz°

08

p-0.8 P - 0 . 6

P,O.4

P-0.2

p=o

F i g . 3. The effect o f correlation (p) and drift (c) on the m e a n recurrence time of emptiness ( m o o ) o f a standard reservoir of size c = 4, fed by a discrete +-1 input (full lines) and

normal input

(dashed lines).

295

We find that for e ( 0.6 very fair approximations obtain, as can be seen in Fig. 2. Fig. 3 shows moo as a function of e and p for c = 4 for both the normal (with figures about 5% accurate) and discrete -+1 input. The discrete reservoir tends to exaggerate the probability of failure, but reflects the general behaviour of the correlated normal input quite satisfactorily.

CONCLUSIONS

(1) The concept of the standard reservoir is useful in that it allows fair comparisons between what are quite disparate input distributions. The results show that storage probabilities are remarkably insensitive to marginal distribution over most of the parameter space. (2) Because of the large sampling variation of the mean time of first return to emptiness, moo, considerable care must be given to interpreting results of simulations, especially as they are so easy to do; to generalize from the results of a few isolated simulations is unwise. (3) From the results of the few computations reported herein, the m e t h o d of numerical solution of the integral equation describing the storage distribution function is more efficient than the difference equation method. Although it is more difficult to set up in the dependent input case, it is an avenue worth pursuing, especially as error bounds on the resulting computations can be easily found. (4) The random walk approximation to the reservoir problem is very easily computed, and yields fair estimates of moo and m0c even though the normal and discrete inputs are so different. We do, however, find that for meaningful results to come from discrete reservoir approximations, the "size" k must be carefully related to the standard capacity of the continuous reservoir. If this is done, then other discrete net input distributions {such as Poisson) should yield better results than the + 1 input, especially if skewness is explicitly taken into account. (5) We find that an increase in serial correlation, p, decreases moo , as expected. (6) We find that an increase in skewness, 7, may decrease moo (for small e and c), but in general increases moo (when c is large and e > 0). This effect is delayed when p is large. (7) We find that a great deal of computing power is required when c, ~ and p are large, even for moderate e, when we wish to ascertain moo. It is felt that the methods of numerical solution of difference and integral equations may n o t be the only methods which are confounded by this difficulty - - i t should be borne in mind when conducting simulation experiments. (8) The numerical results, which are accurate to three significant digits or more, should provide useful checks for various methods of determining reservoir reliability. It appears that a useful contribution could be made by fitting approximating functions to the results reported here, as they represent a large computing investment.

296 REFERENCES Baker, C.T.H., 1977. The Numerical Treatment of Integral Equations. Clarendon, Oxford. Feller, W., 1968. An Introduction to Probability Theory and Its Applications, Vol. 1. Wiley, New York, N.Y., 3rd ed. Gomide, F.L.S., 1975. Range and deficit analysis using Markov chains. Colo. State Univ., Fort Collins, Colo., Hydrol. Pap. 79. Johnson, N.L. and Kotz, S., 1972. Distributions in Statistics: Continuous Multivariate Distributions. Wiley, New York, N.Y. Kartvelishvili, N.A., 1969. Theory of Stochastic Processes in Hydrology and River Runoff Regulation. Israel Program for Scientific Translations, Jerusalem (translated from Russian). Kemeny, J.G. and Snell, J.L., 1960. Finite Markov Chains. van Nostrand, London. Lloyd, E.H., 1963. A probability theory of reservoirs with serially correlated inputs. J. Hydrol., 1: 99--128. Moran, P.A.P., 1959. The Theory of Storage. Methuen, London. Pegram, G.G.S., 1974. Factors affecting the probability distribution of draft from a Lloyd reservoir. Water Resour. Res., 10: 63--66. Pegram, G.G.S., 1975. Recurrence times of draft-patterns from reservoirs. J. Appl. Prob., 12: 647--652. Pegram, G.G.S., 1978. Some simple expressions for the probability of failure of a finite reservoir with Markovian input. Geophys. Res. Lett., 5: 13--15. Sowden, R.R. and Ashford, J.R., 1969. Computation of the bivariate normal integral. Appl. Stat., 18: 169--180.