Journal of Hydrology 2 (1964) 1-10; © North-Holland Publishing Co., Amsterdam N o t t o be reproduced by photoprint or microfilm without written permission from the publisher
PROBABILITY THEORY OF RESERVOIRS WITH
SEASONAL I N P U T E. H. LLOYD and S. ODOOM Imperial College, London Abstract: After outlining the major assumptions in the simple probability theory of reservoirs, the paper discusses the modifications required when the inflows are seasonal, and gives a detailed treatment of a simple 2-season model.
1. Introduction
The sequence of levels of a real reservoir is influenced by the (usually) complicated pattern of inflows, by the withdrawal policy, and by what we may call secondary factors such as evaporation losses, blow-off, seepage, etc. A simple mathematical model that may serve as a first approximation to reality can be constructed by adopting the following drastic assumptions: the so-called secondary factors may be neglected; time is discrete, not continuous; withdrawals are instantaneous, with one withdrawal in each time unit; the input pattern forms a fixed probability distribution, unchanging with time; inflows during non-overlapping time intervals are mutually independent. It follows from these assumptions that the sequence of reservoir levels possesses a serial correlation structure known as a Markov Chain. The operating theory of such reservoirs was initiated by Moran 1), and applications have been studied by (amongst others) Gani 2), Prabhu a), Langbein 4), Lloyd and OdoomS). The restriction to discrete time and instantaneous withdrawals need not, in principle, be a serious one, since the units of time may be made very small. In practice it must be admitted that a fine time-scale would increase both the statistical and the computational complexity of the model. A certain amount of progress has been made with a theory of continuous time and continuous withdrawals (Gani and Prabhu6)). The realism of the simple model can, of course, be improved in various ways. One improvement is to abandon the requirement of independent inputs; for example, the inflows may be supposed to possess a Markovian 1
2
E . H . LLOYD AND S. ODOOM
serial structure 7). A n o t h e r improvement would be achieved by relaxing the requirement that the inflow distribution is unchanging in time. The present paper is concerned with this aspect; we suppose in fact that the inflow pattern is seasonal. 2. Seasonal inputs
Suppose we have a reservoir f r o m which m units o f water are removed in each o f the k seasons into which the year is divided. During season 1 o f year t, a total quantity X1, , o f water enters the reservoir, and, at the end of the season, m units are removed (or, at least, as m a n y units as possible are removed, up to a m a x i m u m o f m). It is convenient to regard this removal as being effectively instantaneous. The level o f water is then observed; call it ZI,,. During the second season o f year t, the inflow is X2,t and the level at the end o f the season, after the withdrawal has been made, is Z2, ,. The process goes on in this way, producing a sequence that m a y be illustrated thus: Inflow during Level at end Year Season season of season t
t + 1
t+ 2
1 2 :
XI,t X2,t i
ZI,t Z2,t :
k
X~,t
Z~,t
1
Xl,t+l :
Z1, t+l :
k
Xk,t+l
Z~, t+l
1 !
Xl, t+~ :
Zl,t+~ i
k
X~, t+~
Z~, t+2
Instead o f assuming a fixed distribution, the same for all inflows, as in the non-seasonal model, we now assume a fixed "first-season" distribution, a fixed "second-season" distribution and so on. In the present treatment we suppose all inputs to be mutually independent. We denote the first-season probabilities, (those o f X 1 , t, X 1 , t + 1, X ! ,t + 2 . . . . ) by P(Xl,r+h=r)=pl(r), r = 0 , I , 2 .... f o r h = 0, 1,2 . . . . Likewise, for the second season, let P(X2,,+ h = r) = p2 (r), and so on.
PROBABILITY THEORY OF RESERVOIRS WITH SEASONAL INPUT
3
The object of the exercise is to obtain information about the distribution of levels. As in the case of inflows, the levels Z1, t, Zl,t+l, Zl,t+2 . . . . will have a common distribution. Let
P(Zl,,+,,=s)=(l(s),
s = O, 1,..., n.
P(Z2,t+h = s) = (2(s),
s = 0, 1..... n,
Similarly, let
and so on. Here the maximum value of s is determined by the size of the reservoir. By well-known arguments (see e.g. Lloyd s)) it is easy to derive a relationship between the (2's and the (l's, of the form
(2 (r) = E b~2) (1 (s), 8
where the coefficients O~.(2) the so-called transition probabilities, are combirs nations of the P2'S. Likewise (a (r) = E b~a, (2 (s) $
where the b(a),s are combinations of the pa's. Continuing this process, we ultimately arrive at S
and
(1 (r) = E b~2)(k(S) • 8
The whole set of these relations can be telescoped together; for example, the first and second give (on eliminating the (2's) (3 (~) = E C,dl (s) s
where Crs ~
~-~ b(3)h(2)
~
rj
~js
•
J
The completely telescoped set gives a relation of the form
8
This is a set of (homogeneous) simultaneous linear equations in the (~'s. The solution gives the distribution of levels at the end of the annual "first seasons". The annual r-th seasons may be similarly dealt with.
4
E. H . L L O Y D A N D S. ODOOM
3. A simple case
I n a numerical study the formation o f the coefficients in the telescoped process would be a standard computer operation. In the present note, however, we propose to illustrate the argument by using a very simple example, namely a 2-season year, with 3-valued input distributions. This drastic simplification makes it possible to obtain explicit expressions for the distribution o f levels in each season. The assumption o f a 3-valued input is a formal device corresponding to the a d o p t i o n o f a coarse grouping o f the actual inputs (which in reality o f course will have a continuous range o f possible values). With this convention we suppose the inputs to have possible values m - a, m, m + a, where m is the withdrawal rate. Let the distribution of inputs in the first and second seasons, respectively, be P(X l=m-a,m,rn+a)=(p0,pa,p2) and P(X 2
= m -- a, m, m + a) = (qo, qx, q2).
The reservoir has size n a + m , and thus has effective capacity na. The levels at the end o f the first and second seasons are denoted by Z1, Z2 respectively, and their distributions by P(Z 1 = ra) = f,,
r = O, 1 . . . . , n
P(Z2 = sa) = gs,
s = O, 1 . . . . . n .
and
The matrix o f transition probabilities f r o m season '2' to season '1' is readily seen (c.f. Langbein4)) to be as follows: ~
=
.
Level at end of 2nd season of year t 0 p0 + pl
1 po
2
0
3
1
p2
pl
po
pz
pl
p0
p2
pl
...
n--1
n
c~
2 3
n--2
po
n-- 1 n
pl
po
p2
pa +p2
PROBABILITY THEORY OF RESERVOIRS WITH SEASONAL INPUT
5
This means that fo = (Po + Pl)go q- Pogl ./'1 = P2go + Pig1 + Pog2 f2 = P2gx + Plgz + Pog3,
etc. The transition matrix from season '1' to season '2' is of the same form, but with q's replacing p's, and with f ' s and g's interchanged, whence go = (qo + q l ) f o + q o f l gl = qzfo + q l f l + qof2 g2 = q z f l + q l f 2 + qof3,
etc. 'Telescoping' these two systems is equivalent to eliminating either t h e f ' s or the g's; if the latter, we obtain the following set of homogeneous linear equations in t h e f ' s : {(Po + Pl)(qo + ql) + Poq2 - 1}fo + {(Po + Pl)qo + Poql} f l + Poqof2 = 0,
(3.1)
{P2 (qo + ql) + P l q 2 } f o + {P2qo + P l q l + Poq2 -- 1}f, + {Pxqo + Poql}f2 + Poqofa = 0,
(3.2)
P2q2f, + {P2ql + P l q 2 } f , + l + {P2qo + P l q l + Poq2 -- 1} f,+2 + {Plqo + P o q l } f , + 3 + Poqof~+4 = O,
r = 0, 1.... , n -- 4,
(3.3)
P 2 q 2 f . - a + {P2ql + Plq2} f n - 2 + {P2qo + Pxql + Poq2 -- 1} f . - t + {Plqo + Po(ql + q2)} f . = 0, (3.4)
P2q2f.-2 + {P2q2 + q2(Pl + P2)} f . - 1 + {P2qo + (Pl + P2)(ql + q2) -- 1}f, = 0.
(3.5)
It is helpful to regard this set as being made up of the difference equation (3.3), with side conditions (3.1), (3.2), (3.4), and (3.5). Equation (3.3) is a linear difference equation with constant coefficients (the discrete analogue of a linear differential equation with constant coefficients), of the fourth order. It belongs to a class having a general solution of the form fr = A~).~ + A22~ + A32~ + A # ~ ,
r = 0, 1,..., n,
where the A's are constants to be determined by using the side conditions, and the ).'s are the four roots of the quartic P2q2 + ( P 2 q l .-I- plq2)). -I- (P2qo + Plq~ + P o q 2 -- 1) •2
+ (Plqo + Poqx) )-3 + Poqo 2-4 = 0,
6
E.H. LLOYD AND S. ODOOM
provided these roots are distinct. It is easily verified that 2 = 1 is a root, so that the solution reduces to f~ = At2~ + A22[ + A32~ + A4
(3.6)
where the 2's are now the roots of the cubic Poqo 23 + {Po (qo + ql) + Plqo} 2 2 - {Plq2 + P2 (ql + q2)} 2 - p 2 q 2 =0. (3.7) On substituting the solution (3.6) in the 'side condition equations' it will be found that A4 = 0, while A 2 and A a may be expressed in terms of A t by the equations
2223 -- PO A1 +
2321
2~ 21-qo/At + 2"2
Po/
Po/
2~-qo/A2+2; 22
go/
Finally A 1 is fixed by the condition ~f,=l
r=O
where f , is given by (3.6). With these formulae it would be easy to tabulate the probabilities f , associated with all possible first-season levels for any choice of first-season and second-season inflow probabilities; and, since the g's are known linear functions of the f ' s , the second-season values may then also be simply computed. 4. Numerical studies
It is not the object of the present note to embark on large-scale tabulation, and we content ourselves with numerical treatment of a special case of the situation discussed in section 3. The special case in question is that which arises when q0 = P2
and
q2 ----PO,
(Po ~ P 2 ) ,
that is, when the input distributions in the two seasons are not individually symmetrical and are mirror-images of one another, in the sense that t(Xl=m-a,m,m+a)=(po, t(X2=m-a,m,m+a)=(p2,
pl, p2) pl, po).
Sincepo + P l q-P2 = 1, the solution here will involve only two parameters. It turns out for this case that 2 = 1 is a double root of the quartic for 2 (i.e.
PROBABILITY THEORY OF RESERVOIRS WITH SEASONAL INPUT
7
the cubic (3.7) that remains after extracting the root 2 = 1 is itself satisfied by 2 = 1). When there are coincident roots the general solution of a difference equation takes a different form, (see, for example, Levy and Lessman 9)) and in the present case (3.6) is replaced by f , = Ai2 ~ + A22r2 + A 3 + A4r ,
r = O, 1 ..... n,
where 21 and 22 satisfy the quadratic that remains after extracting the second unit root. This is
22 +(POP 5 + PxP2 + 2 ) 2 - 1 = 0 . \
PoP2
The side equations lead to A4 = 0 and to the following equations for A t, A2 and A3 : ( 2 2 - P2)At + ( 2 1 - P 2 ~ A 2 + ( 1 - P 2 ~ A a = Po/ Po/ Pool
0
2 ] ( 21-p~°~A1p2/ + 2 " 2 ( 2 2 - ~ 2 ) A 2 + ( I - P ° ~ A a = O p 2 / whence we find f, = k {(P2 -- Po21)(2~ +1 -- 1)2] -- (P2 -- po22)(2] +1 -- 1)2~ + ct} r = 0 , 1 ..... n where _ po + p
(2y1 _ 2 ] + , )
Po - P2 and k is the normalising constant defined by ~f~=
1.
0
The second-season distribution is obtained from this by merely interchanging Po and P2. The roots 21, 22 are not affected by this interchange; the constant e changes sign. In order to display the dependence of the f , on the parameters of the problem we introduce the average inputs 121 and 122 for the two seasons respectively: I~ = E ( x l ) = m + (P2 - Po)a P2 = E(x2) = m - (P2 - Po) a and their common variance a 2 = var(xl) = var(x2) ={P2+Po-(P2-Po)
2}a2-
8
E . H . LLOYD AND S. ODOOM
We propose to use Y=
~ (P2 - Po)
as a dimensionless measure of the difference between the inflow distributions of the two seasons. When y = 0 the two seasons have identical input distributions, so that we are back to a single-season problem; further, the common inflow distribution is then symmetrical, and in that case the reservoir levels are known to be uniformly distributed s). Increasing values of y correspond to increasing skewness of the inflow distribution in each season, the skewnesses being, of course, in opposite directions in the two seasons. To tabulate the distribution of reservoir levels for various values of y, it is necessary to fix the value of a/cr. Since a three-valued input can only be a convenient formal device for describing a real input having a continuum of values, the choice of value for a corresponds to the choice of spacing of the three ordinates which are being used to approximate the continuous real distribution of inflows. Following Lloyd s) we choose a = 1.5 or. We then obtain the following table relating values of y, Po, P2 and 21. [-The form of the quadratic equations for 21, 22 shows that 2122 = 1. We henceforth use 21 to denote the root such that 1211 > 1 ; then 22 = 1/21]. TABLE 1 y
po
p~
21
0.1
.191
.258
- - 6.88
0.2
.164
.298
- - 6.93
0.3
.142
.342
- - 7.00
0.4
.124
.391
- - 7.01
0.5
.111
.444
- - 6.84
0.6
.102
.502
- - 6.53
0.7
.098
.564
- - 5.87
0.8
.098
.631
- - 4.99
0.9
.102
.702
- - 3.95
1.0
.111
.778
- - 2.77
F o r c o r r e s p o n d i n g n e g a t i v e v a l u e s o f y , interc h a n g e po a n d p z ; 2 r e m a i n s u n c h a n g e d .
The significance of the sign of y is that a positive y corresponds to a "wet" Season 1 and a " d r y " Season 2; for a negative y the positions are interchanged; larger numerical values o f y correspond to greater differences between the seasons. The next table shows the probabilities of emptiness, for various values of y:
9
PROBABILITY THEORY OF RESERVOIRS WITH SEASONAL INPUT TABLE 2 Probability of emptiness at end of season
y
n=10
-- 1.0 --0.8 --0.6 --0.4 --0.2 0 0.2 0.4 0.6 0.8 1.0
.120 .108 .104 .102 .097 .091 .083 .074 .061 .049 .036
Season 1
Season 2
Size of Reservoir
Size of Reservoir
n
15 .082 .074 .071 .070 .066 .063 .057 .051 .042 .037 .025
n~20
n=
.062 .056 .054 .053 .051 .048 .043 .038 .032 .025 .019
10
.036 .049 .061 .074 .083 .091 .097 .102 .104 .108 .120
n=
15
n =20
.025 .037 .042 .051 .057 .063 .066 .070 .071 .074 .082
'019 .025 .032 .038 .043 .048 .051 .053 .054 .056 .062
O f p o s s i b l y m o r e p r a c t i c a l i n t e r e s t is t h e p r o b a b i l i t y o f a t l e a s t o n e e m p t i n e s s d u r i n g t h e y e a r . T h i s is t h e s u m o f t h e t w o c o r r e s p o n d i n g p r o b a b i l i t i e s i n t h e a b o v e t a b l e . T h e v a l u e s a r e as s h o w n i n T a b l e 3: TABLE 3 Probability of at least one emptiness during the year y
n = 10
n = 15
n =20
1.0 0.8 -0.6 -0.4 -- 0.2 0 0.2 0.4 0.6 0.8 1.0
.156 .157 .165 .176 .180 .182 .180 .•76 .165 .157 .156
.107
.081 .081 .086
--
--
.111
.113 .121 .•23 .126 .123 .121 .113 .111 .107
.091
.094 .096 .094 .091 .086 .081 .081
A n i n t e r e s t i n g f e a t u r e o f t h i s t a b l e is t h a t , w h i l e t h e p r o b a b i l i t y o f a t l e a s t o n e e m p t i n e s s is s o m e w h a t g r e a t e r w h e n t h e r e is n o d i f f e r e n c e b e t w e e n t h e two seasons than in the case of true seasonality, the variation in this probab i l i t y w i t h y is c o m p a r a t i v e l y slight. I n o t h e r w o r d s , f o r t h e 2 - s e a s o n m o d e l
10
E. H. LLOYD AND S. ODOOM
investigated here, the probability of getting at least one emptiness does n o t seem to depend very strongly o n the degree of difference between the seasons. Whether this is a peculiarity of the particular case examined is n o t clear; further studies of this question will be u n d e r t a k e n .
References 1) P. A. P. Moran, The Theory of Storage (Methuen and Co. Ltd., London, 1959) 2) J. Gani, Problems in the Probability Theory of Storage Systems, J. R. Statist. Soc. B. 19 (1957) 181 3) N. U. Prabhu, Q. J. Maths. 9 (1958) 4) W. B. Longbein, Reservoir Storage - General Solution of a Queue Model. Geological Survey Research (1961) art. 298 5) E. H. Lloyd and S. Odoom, A Note on the Solution of Dam Equations, J.R.S.S. (B) (1963): in the press. 6) J. Gani and N. U. Prabhu, A Storage Model with ContinuousInfinitelyDivisible Inputs, Proc. Camb. Phil. Soc. 59 (1963) 417 7) E. H. Lloyd, Reservoirs with Serially Correlated Inflows, Technometrics 5 (1963) 85 8) E. H. Lloyd, A Probability Theory of Reservoirs with Serially Correlated Inputs, J. of Hydrology 1 (1963) 99 9) H. Levy and F. Lessman, Finite Difference Equations (Pitman, London, 1959)