# The Probability of Success of Certain Storage Schemes

A.-A. Anis

► To cite this version:

A.-A. Anis. The Probability of Success of Certain Storage Schemes. Annales de l'ISUP, 1965, XIV, pp.407-422. hal-04084895

HAL Id: hal-04084895

https://hal.science/hal-04084895

Submitted on 28 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# THE PROBABILITY OF SUCCESS OF CERTAIN STORAGE SCHEMES

## A.-A. ANIS

I. Most of the papers on the theory of storage, as developed by Moran ans his school, are centred around the question of the distribution of the storage content $t_n$, which is defined at discrete times $n = 0, 1, 2, ...$ by the recurrence relations

$$t_{n+1} = \min(k, t_n + x_n) - \min(m, t_n + x_n), \quad 0 < m < k \quad (1.1)$$
where

(i) $x_n$ denotes the dam input during the n-th year, and it is assumed that $x_0, x_1, ...$ are mutually independent and identically distributed variates with a finite mean $E(x)$.

(ii) $k$ is the capacity of the dam, and thus the term $\min(k, t_n + x_n)$ represents a possible overflow which occurs if and only if $t_n + x_n > k$.

(iii) the term $\min(m, t_n + x_n)$ indicates a release policy according to which, at time $(n+1)$, an amount $m$ of water is released, unless the dam contains less than $m$, in which case the entire available amount is released.

Instead of studying the distribution of $t_n$, we propose to discuss directly methods of obtaining the probability of the "success" of any defined release rule, where success is defined as availability of water whenever it is required for release.

We start by making the two assumptions :

(i) the $x_n$'s are independent variates with identical probability distributions and the mean $E(x)$ is finite.

(ii) $y_1, y_2, ... y_n$ are the random outflows, and these $y_n$'s are also independent variates with identical probability distributions and the mean $E(y)$ is finite.

The choice of the outflow as random generalises Moran's model, and does not of course exclude the possibility of a regular demand $m$. We will work at the end of the paper an example in provisioning theory to show how our general equations work.
We shall also make the assumption that the scheme starts with an initial amount of water $t_0$ where $0 \le t_0 \le k$

Now during the first year the dam receives the random inflow $x_1$ and releases the random outflow $y_1$; thus ending with a quantity $t_1$ in the dam where

$$t_1 = \min(k, t_0 + x_1 - y_1)$$

During the second year, the dam receives the random inflow $x_2$ and releases the random outflow $y_2$, thus ending with $t_2$ where

$$t_2 = \min(k, t_1 + x_2 - y_2)$$

and so on

In general

$$t_n = \min(k, t_{n-1} + x_n - y_n), \quad n \ge 1$$

We stipulate that the scheme will "succeed" in the time-period $n$ if and only if $t_1, t_2, \dots, t_n$ are all non-negative. In fact this will be our definition of success throughout the paper.

Now let us define a new set of variates $u_1, u_2, \dots, u_n$ where,

$$u_r = y_r - x_r$$

$$r = 1, 2, \dots, n$$

and $G(u)$ is the distribution function of the u's

Also let $P_n(t)$ denotes the probability that starting with an initial quantity $t$ in the dam, the scheme will succeed during the following $n$ years. If $x$ is the quantity received during the first year, $y$ the quantity demanded for release, the scheme would fail if $x + t - y < 0$

Thus for success

$$x + t - y \ge 0$$

or

$$u \le t$$

leaving in store $t - u$ if this quantity is less than $k$. If it is greater than $k$, then it will leave $k$ itself and the excess is wasted.

Hence if

(a) $t - u < 0$, the scheme fails

(b) $0 \le t - u < k$, the scheme succeeds with a probability
$P_{n-1}(t - u)$

(c) $t - u \ge k$, the scheme succeeds with a probability
$P_{n-1}(k)$

The relation between $p_n, p_{n-1}$ is given by the equation

$$P_n(t) = \int_{0 \le t-u \le k} P_{n-1}(t-u) \, dG(u) + \int_{t-u \ge k} P_{n-1}(k) \, dG(u) \quad (1.2)$$
Where the integration sign denotes Stieltjes integral, so that we can combine discrete and continous cases in a single equation.

Equation (1.2) is a recurrence relation which, theoretically speaking, could be used to obtain $P_n(t)$ successively starting with the initial condition $P_0(t) = 1$. In few cases it is possible to obtain a general formula for $P_n(t)$, but in the most important cases this is not the situation.

It is worth noting, on passing, that $P_n(t)$ depends only on the distribution of the difference between the outflow and the inflow, and not on their individual distributions. This property seems intuitive, and corresponds to the fact, in the theory of queues with a single server, that the distribution of waiting times depends on the distribution of the differences between serving times and arrival times and not on their individual distributions.

## II. Let us take, as an example of equation (1.2), the case of a laplace input distribution

$$dF(x) = \frac{1}{2} e^{-|x-\mu|} dx \quad (2.1)$$

where there is a regular demand $\mu$ every year.

In this case (1.2) becomes

$$P_n(t) = \frac{1}{2} \int_0^k P_{n-1}(w) e^{-|w-t|} dw + \frac{1}{2} P_{n-1}(k) \int_k^\infty e^{-|w-t|} dw$$

ore more precisely

$$P_n(t) = \frac{1}{2}e^{-t} \int_0^t P_{n-1}(w) e^w dw + \frac{1}{2}e^t \int_k^t P_{n-1}(w) e^{-w} dw \quad (2.2) \\ + \frac{1}{2} P_{n-1}(k) e^{t-k}$$

We define the probability generating function of $P_n$'s as

$$\psi(t) = \sum_{r=0}^{\infty} \lambda^r P_r(t) , \quad |\lambda| < 1 \quad (2.3)$$

This is a uniformly convergent series, and if we multiply (2.2) by $\lambda^n$ and sum both sides with respect to $n$ from $n=1$ to $\infty$ we get

$$\psi(t) = 1 + \frac{\lambda}{2} e^{-t} \int_{0}^{t} \psi(w) e^{w} dw + \frac{\lambda}{2} e^{t} \int_{t}^{k} \psi(w) e^{-w} dw + \frac{\lambda}{2} \psi(k) e^{t-k} (2.4)$$

On differentiating (2.4) twice we easily get the second order differential equation

$$\psi'(t) = (1 - \lambda) \psi(t) - 1 \quad (2.5)$$
together with the boundary conditions

$$ \left. \begin{array}{l} \psi(0) - \psi'(0) = 1 \\ \psi'(k) + (1 - \lambda) \psi(k) = -1 \end{array} \right\} \qquad (2.6) $$

The general solution of (2.5) would be given by

$$ \psi(t) = ae^{t\sqrt{1-\lambda}} + be^{-t\sqrt{1-\lambda}} + \frac{1}{1-\lambda} \qquad (2.7) $$

where

$$ \begin{aligned} a &= -\lambda(1 - \sqrt{1-\lambda}) / (1-\lambda) \left\{ (1-\sqrt{1-\lambda})^2 + (1+\sqrt{1-\lambda})^2 e^{2\kappa\sqrt{1-\lambda}} \right\} \\ b &= -\lambda(1 + \sqrt{1-\lambda}) / (1-\lambda) \left\{ (1-\sqrt{1-\lambda})^2 + (1+\sqrt{1-\lambda})^2 e^{2\kappa\sqrt{1-\lambda}} \right\} \end{aligned} \qquad (2.8) $$

Thus we finally get

$$ \psi(t) = 1/(1-\lambda) - (\lambda/1-\lambda) \left\{ (1-\sqrt{1-\lambda}) e^{t\sqrt{1-\lambda}} + (1+\sqrt{1-\lambda}) e^{(2k-t)\sqrt{1-\lambda}} \right\} \\ \left\{ (1-\sqrt{1-\lambda})^2 + (1+\sqrt{1-\lambda})^2 e^{2k\sqrt{1-\lambda}} \right\} \qquad (2.9) $$

To obtain explicit expressions for $P_n(t)$, we propose to divide both the numerator and the denominator of the second term of the R.H.S. of (2.9) by $(1 + \sqrt{1-\lambda})^2 e^{2\kappa\sqrt{1-\lambda}}$, thus getting $\psi(t)$ in the form

$$ \psi(t) = \frac{1}{1-\lambda} - \frac{\lambda}{(1-\lambda)(1+\sqrt{1-\lambda})} \cdot \frac{[\alpha e^{-(2k-t)\sqrt{1-\lambda}} + e^{-t\sqrt{1-\lambda}}]}{[\alpha^2 e^{-2\kappa\sqrt{1-\lambda}} + 1]} $$

where

$$ \alpha = \frac{1 - \sqrt{1 - \lambda}}{1 + \sqrt{1 - \lambda}} $$

On expanding $(1 + \alpha^2 e^{-2\kappa\sqrt{1-\lambda}})^{-1}$ we finally obtain

$$ \psi(t) = 1/(1-\lambda) - (1/(1-\lambda) - 1/\sqrt{1-\lambda}) L(t) \qquad (2.10) $$

where

$$ L(t) = e^{-t\sqrt{1-\lambda}} - \alpha^2 e^{-(2k+t)\sqrt{1-\lambda}} + \alpha^4 e^{-(4k+t)\sqrt{1-\lambda}} - \dots \\ + \alpha e^{-(2k-t)\sqrt{1-\lambda}} - \alpha^3 e^{-(4k-t)\sqrt{1-\lambda}} + \dots \qquad (2.11) $$

and thus to get $P_n(t)$ we need to know $L^{(0)}$, $L^{(1)}$, ..., $L^{(n)}$, where $L^{(r)}$ is the rth derivative of $L(t)$ with respect to $\lambda$ at the origin. We notice that $\alpha = 0$ at $\lambda = 0$, and hence in getting $L^{(r)}$ we need not consider the terms of (2.11) which contain $\alpha^{r+1}$ or higher powers. If $e^{-t\sqrt{1-\lambda}}$ is denoted as $\beta$, then the derivatives of $\alpha$ and $\beta$ at the origin can be computed from the obvious relations.
$$ 
\begin{aligned}
\lambda \alpha' &= (1 - \lambda)^{-1/2} \\
\beta' &= \frac{t}{2} (1 - \lambda)^{-1/2} \beta
\end{aligned}
\qquad (2.12)
$$

and on using Libnitz formulae we get

$$ 
\begin{aligned}
(n-1) \alpha^{(n)} &= \sum_{r=1}^{n-1} \left( \frac{2r!}{2^{2r} r!} \right) \binom{n}{r} \alpha^{(n-r)} \\
\beta^{(n+1)} &= \frac{t}{2} \sum_{r=0}^{n} \left( \frac{2r!}{2^{2r} r!} \right) \binom{n}{r} \beta^{(n-r)}
\end{aligned}
\qquad (2.13)
$$

Thus we get

$$ L^{(0)} = e^{-t} $$

$$ L^{(1)} = \frac{t}{2} e^{-t} + \frac{1}{4} e^{-(4k-t)} $$

$$ L^{(2)} = \frac{1}{4} (t^2 + t) e^{-t} + \left\{ \frac{(2k - t) + 1}{4} \right\} e^{-(2k-t)} - \frac{1}{8} e^{-(2k+t)} $$

$$ L^{(3)} = \frac{1}{8} (t^3 + 3t^2 + 3t) e^{-t} $$

$$ 
\begin{aligned}
&+ \frac{1}{32} \left\{ 6(2k-t)^2 + 18(2k-t) + 15 \right\} e^{-(2k-t)} && (2.14) \\
&- \frac{1}{16} \left\{ 3(2k+t) + 6 \right\} e^{-(2k+t)} - \frac{3}{32} e^{-(2k-t)}
\end{aligned}
$$

$$ 
\begin{aligned}
L^{(4)} &= \frac{1}{16} (t^4 + 6t^3 + 15t^2 + 15t) e^{-t} + \frac{1}{16} \left\{ \begin{aligned} &2(2k-t)^3 + 16(2k-t)^2 \\ &+ 31(2k-t) + 21 \end{aligned} \right\} e^{-(2k-t)} \\
&- \frac{1}{16} \left\{ \begin{aligned} &5(2k+t)^2 + 17(2k+t) \\ &+ 21 \end{aligned} \right\} e^{-(2k+t)} \\
&- \frac{1}{16} \left\{ \begin{aligned} &3(4k-t) + 9 \end{aligned} \right\} e^{-(4k-t)} + \frac{3}{32} e^{-(4k+t)}
\end{aligned}
$$

The relations between the P's and the L's are given by

$$ P_1(t) = 1 - \frac{1}{2} L^{(0)} $$

$$ P_2(t) = 1 - \frac{5}{8} L^{(0)} - \frac{1}{2} L^{(1)} \qquad (2.15) $$

$$ P_3(t) = 1 - \frac{11}{16} L^{(0)} - \frac{5}{8} L^{(1)} - \frac{1}{4} L^{(2)} $$

$$ P_4(t) = 1 - \frac{93}{128} L^{(0)} - \frac{11}{16} L^{(1)} - \frac{5}{16} L^{(2)} - \frac{1}{12} L^{(3)} $$

$$ P_5(t) = 1 - \frac{193}{256} L^{(0)} - \frac{93}{128} L^{(1)} - \frac{11}{32} L^{(2)} - \frac{5}{48} L^{(3)} - \frac{1}{48} L^{(4)} $$

...... etc.
In particular if we start with a full dam, we get

$$P_1(k) = 1 - \frac{1}{2} e^{-k}$$

$$P_2(k) = 1 - \frac{1}{4} (k + 3) e^{-k}$$

$$P_3(k) = 1 - \frac{1}{32} (2k^2 + 14k + 29) e^{-k} + \frac{1}{32} e^{-3k}$$

$$P_4(k) = 1 - \frac{1}{192} (2k^3 + 24k^2 + 111k + 195) e^{-k} + \frac{1}{128} (3k + 10) e^{-3k} \quad (2.16)$$

$$P_5(k) = 1 - \frac{1}{768} (k^4 + 18k^3 + 142k^2 + 532k + 843) e^{-k} + \frac{1}{152} (30k^2 + 55k + 67) e^{-3k} - \frac{1}{512} e^{-5k} \quad \dots\dots\dots \text{etc.}$$

It may be noticed, on passing, that when $k = \infty$, $\psi(t)$ becomes

$$\psi(t) = \frac{1}{1 - \lambda} - \left[ \frac{1}{1 - \lambda} - \frac{1}{\sqrt{1 - \lambda}} \right] e^{-t \sqrt{1-\lambda}} \quad (2.17)$$

$$\text{and thus } P_n(o) = (2n)! / 2^{2n} (n!)^2, \quad (2.18)$$

which is the coefficient of $\lambda^n$ in the expansion of $(1 - \lambda)^{-1/2}$.

This is the same result if the input distribution was normal instead of Laplace. This may be seen from the fact that $P_n(t)$ defined by the recurrence relation

$$P_n(t) = \int_0^\infty P_{n-1}(w) \frac{e^{\frac{1}{2}(w-t)^2}}{\sqrt{2\pi}} dw, \quad (P_0 = 1) \quad (2.19)$$

obeys the identity

$$\sum_{r=0}^{n} P_{n-r}(o) P_r(o) = 1 \quad (2.20)$$

[See Anis (1)]

And on applying the method of generating functions on (2.20), we get (2.18).

It should be mentioned that the special case of the Laplace input, when $k = \infty$, was investigated separately by Simaika [2], but his approach is of a different nature.
### III. Let us now study the case of a discrete input, when k is finite. We will assume, for the sake of simplicity that k, t, and the x's are all expressed as integers, and that the release is regular and is taken as unity.

Equation (1.2) thus takes the form

$$P_n(t) = \sum_{x+t-l=0}^{k-1} P_{n-1}(x+t-1) p_x + P_{n-1}(k) \sum_{x+t-1=k}^{\infty} p_x \quad (3.1)$$

where $p_x$ is the probability distribution of the input.

This equation may be put in matrix form. For let $\underline{P}_n$ denotes the vector given by

$$\underline{P}'_n = \{P_n(0), P_n(1), P_n(2), \dots, P_n(k)\}$$

Then equation (3.1) gives us

$$\underline{P}_n = \underline{\Delta} \underline{P}_{n-1} \quad (3.2)$$

where $\underline{\Delta}$ is the $(k + 1 \times k + 1)$ matrix

$$\underline{\Delta} = \begin{bmatrix} p_1 & p_2 & p_3 & \cdots & p_k & q_k \\ p_0 & p_1 & p_2 & \cdots & p_{k-1} & q_{k-1} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ p_0 & p_0 & p_1 & \cdots & p_{k-2} & q_{k-2} \\ \multicolumn{6}{c}{\text{...........}} \\ p_0 & p_0 & p_0 & \cdots & p_0 & q_0 \end{bmatrix} \quad (3.3)$$

and $q_r$ is given by

$$q_r + \sum_{l=0}^{r} p_l = 1 \quad r = 0, 1, 2, \dots \quad (3.4)$$

Formally the solution of (3.2) is given by

$$\underline{P}_n = \underline{\Delta}^n \underline{1} \quad (3.5)$$

where $\underline{1}$ is the unit vector

Further if the characteristic roots of $\underline{\Delta}$ are distinct or if the degeneracy of the characteristic matrix of $\underline{\Delta}$ is equal to the number of equal roots, then $\underline{\Delta}$ could be put in the diagonal form

$$\underline{\Delta} = \underline{H} \underline{\Delta} \underline{H}^{-1} \quad (3.6)$$

where $\underline{H}$ is the matrix of the characteristic vectors and $\underline{\Delta}$ is the diagonal matrix of the characteristic roots.

In this case

$$\underline{P} = \underline{H} \underline{\Delta}^n \underline{H}^{-1} \underline{1} \quad (3.7)$$
In some cases, the determination of the vector $\underline{P_n}$ is possible even without using this approach.

For example, in the binomial case $p_0 = p, p_1 = q, (p+q=1)$, the matrix $\underline{\Delta}$ takes the form

$$ \underline{\Delta} = \begin{bmatrix} q & o & o & \dots & o & o \\ p & q & o & \dots & o & o \\ o & p & q & \dots & o & o \\ \multicolumn{6}{c}{\dots\dots\dots\dots\dots\dots\dots} \\ o & o & o & p & q & \dots \end{bmatrix} \quad (3.8) $$

and direct calculation of $\underline{\Delta}^n$ would give us the binomial coefficients in the lower half of the matrix. For example, when $k = 3$, $\underline{\Delta}^n$ becomes

$$ \underline{\Delta}^n = \begin{bmatrix} q^n & o & o & o \\ (\binom{n}{1})p q^{n-1} & q^n & o & o \\ (\binom{n}{2})p^2 q^{n-2} & (\binom{n}{1})p q^{n-1} & q^n & o \\ (\binom{n}{3})p^3 q^{n-3} & (\binom{n}{2})p^2 q^{n-2} & (\binom{n}{1})p q^{n-1} & q^n \end{bmatrix} $$

This leads us easily to the formula

$$ P_n(t) = \sum_{r=0}^{t} \binom{n}{r} p^r q^{n-r} \quad (3.9) $$

which may be easily proved by induction on using (3.1).

The case of the geomatric input distribution could be also obtained by methods similar to those of Weesakul [3].

#### IV. Now let us consider the case of an infinite dam, i.e $k = \infty$. In this case equation (1.2) reduces to

$$ P_n(t) = \int_{u \le t} P_{n-1}(t-u) \, dG(u) \quad (4.1) $$

This is exactly the equation obtained by Lindley [4] on the distribution of the waiting time of the $(n+1)$ customer in the theory of queues. This formulation shows that from the point of view of mathematical formalism, the demand here corresponds to the service time in the queue's theory, and the supply corresponds to the arrival times.
This aspect of the theory of storage, as developed here, could be also related to the theory of random walk with a single absorbing barrier. In the case of an infinite dam, the scheme will succeed if

$$u_1 \le t, u_1 + u_2 \le t, \dots, u_1 + u_2 + \dots u_n \le t$$

Hence

$$P_n(t) = P_r \{ u \le t, u_1 + u_2 \le t, \dots, u_1 + u_2 + \dots u_n \le t \}$$

$$= P_r \{ U_n \le t \}$$

where $U_n$ is the maximum of the sequence of partial sums of the u's. Such derivation shows that $P_n(t)$ is a distribution function.

Now consider a particle performing random walk on the x-axis with an absorbing barrier situated at a position $x = t$ ($t \ge 0$). The particle starts its motion at the origin, moving in jumps of amount $u_r$ successively, until the barrier is crossed. $P_n(t)$ then is the probability that after $n$ steps the walk is still in progress, and the particle did not cross the barrier, using of course the convention that the walk which arrives at $t$ will not have crossed the barrier.

Once relation (4.2) is established, we can proceed to use a theorem due to Lindley [4] namely:

The necessary and sufficient condition that $P_n(t)$ has a non-degenerate limit as $n \to \infty$ is that $E(u) < 0$, or $u = 0$ certainly.

The case of $u=0$ certainly will be excluded on account of its triviality, and so we will be assuming, from now onwards, that $E(u)$ is negative unless otherwise stated.

Now suppose that the release demand is constant (m say) and that the mean supply is $\mu$. We have

$$E(u) = m - \mu \quad (4.3)$$

If $E(u) \ge 0$, then $P_n(t) \to 0$, as $n \to \infty$

i.e. if the constant demand is greater or equal to the mean supply, the scheme will fail in the long run with probability one, whatever we initially store. This is only to be expected, in view of the severe conditions we are imposing. Only in the case where the constant demand is less than the mean supply, does the question of limiting probability arise.

Now consider the case $E(u) < 0$,
when $P_n(t) \rightarrow P(t)$, equation (4.1) becomes

$$P(t) = \int_{u \le t} P(t-u) \, dG(u) \qquad (4.4)$$

or

$$P(t) = \int_{w \ge t} P(w) \, dG(t - w) \qquad (4.5)$$

If $G(u)$ is absolutely continuous, so the derivative $g(u)$ of $G(u)$ exists, then (4.5) may be written

$$P(t) = \int_{0}^{\infty} P(w)g(t - w) \, dw \qquad (4.6)$$

This is an integral equation of the Wiener and Hopf type, and provided certain conditions on $g(w)$ are satisfied, it may be solved by the generalisation given by Smithies [5]. However in the cases discussed here, it is easier to use other methods.

As an example, let us take the case of the negative exponential input. This case is interesting, since it is known that in many cases, the use of the negative exponential gives good fit to rainfall data.

It is convenient to take the mean input to be equal to unity, and thus the constant demand $m$ is less than one.

Equation (4.4) gives us

$$P(t) = \int_{x=m-t}^{\infty} P(x+t-m) \, dF(x) , \quad 0 < t < m \qquad (4.7)$$

and

$$P(t) = \int_{x=0}^{\infty} P(x + t - m) \, dF(x) , \quad t > m \qquad (4.8)$$

where

$$\begin{aligned} dF(x) &= e^{-x} dx, \quad x > 0 \\ &= 0, \quad x < 0 \end{aligned} \qquad (4.9)$$

It may be observed that equation (4.7) gives

$$P(t) = P(t) , \quad 0 < t < m \qquad (4.10)$$

giving a solution

$$P(t) = P(0) e^t \qquad (4.11)$$

Also from equation (4.8) one could verify, on using integration by parts, that

$$P(t) = P(t) - P(t-m) , \quad t > m \qquad (4.12)$$
We start with the assumption that $P(t) = 0$, for $t < 0$, jumps from zero to $P(o)$ for $t = 0$, but varies continuously for all values of $t$.

The solution of (4.12) could be obtained by successive substitutions using (4.11)

For example for $m < t < 2m$

$$P(t) = P(t) - P(o) e^{t-m}$$

and on solving this differential equation, we get

$$P(t) = P(o) \{ e^t + (m - t) e^{t-m} \}$$

Proceeding in this way, we finally obtain that when

$$n m < t < (n + 1) m$$

$$P(t) = P(o) \sum_{r=0}^{n} \frac{(rm - t)^r}{r!} e^{t-rm} \quad (4.13)$$

To obtain $P(o)$, we appeal to Jensen's [6] theorem which shows that the sum on the R.H.S of (4.13) tends to $\frac{1}{1-m}$ as $n \to \infty$. Further we know that $P(t) \to 1$ as $n \to \infty$

$$\text{Hence} \qquad P(o) = (1 - m)$$

and thus

$$P(t) = (1 - m) \sum_{r=0}^{n} \frac{(rm - t)^r}{r!} e^{t-rm} \quad (4.14)$$

This result was obtained by Erlang as the distribution of waiting times in the theory of telephone calls with a single operator. The function $P(t)$ was tabulated by Erlang [6].

V. Now suppose that the demand is constant, and for convenience we take it unity in the supply and demand scale. Suppose also that the input is random with a distribution discrete or continuous $F(x)$ say ($x \ge 0$).

$$\text{Thus} \qquad u = 1 - x$$

with the corresponding assumption $E(x) > 1$

$$\text{Since} \qquad G(u) = F(1 - u)$$

equations (4.4) becomes

$$P(t) = \int_{u \le t} P(t - u) dF(1 - u) \quad (5.1)$$

or replacing $t$ by $t + 1$, we get
$$ 
\begin{aligned}
P(t + 1) &= \int_{1-u \ge t} P(t + 1 - u) \, dF(1 - u) \\
&= \int_{w \ge t} P(t + w) \, dF(w)
\end{aligned}
$$

If $t \ge 0$, we finally get

$$ P(t + 1) = \int_{0}^{\infty} P(t + w) \, dF(w) \quad (5.2) $$

Let us deal with the case of discrete inputs, and where the initial store $t$ is also expressed in units.

If we try a solution of the form

$$ P(t) = Z^t, \quad |z| < 1 \quad (5.3) $$

and substitute in (5.2) we get the equation

$$ \Phi(z) = z \quad (5.4) $$

where $\Phi(z)$ is the probability generating function of the input distribution.

Each root of equation (5.4) which satisfies the condition

$$ |z| < 1, $$

gives rise to a solution of (5.2). It is obvious that $Z = 1$ is a root of (5.4). Further $P(t)$ must tend to one as $t$ tends to infinity. This suggests therefore that we should take

$$ P(t) = 1 - \sum_{i=0}^{N} c_i Z_i^t, \quad |z_i| < 1 \quad (5.5) $$

as our solution. The $c_i$s are arbitrary constants, and $N$ is the number of roots of (5.4) which lie inside the unit circle.

**THEOREM.** - If the modulus of the characteristic function $\psi(\theta)$ of the input distribution is less than unity in the range $-\Pi < \theta < \Pi$ (except of course at $\theta = 0$)

then

$$ P(t) = 1 - \alpha^{t+1} \quad (5.6) $$

where $\alpha$ is the unique root of (5.4) which lies inside the unit circle and is real.

We start first by observing that the condition concerning $\psi(\theta)$ is satisfied in familiar discrete distributions. (Bionomial, geometrical... etc).
To prove this theorem, we appeal directly to Rouché theorem in the theory of complex variable. This theorem states :

If $f(z)$ and $g(z)$ are two regular functions within and on a closed contour $C$ on which $f(z)$ does not vanish ; and if

$$|g(z)| < |f(z)|$$

on $C$, then $f(z) + g(z)$ and $f(z)$ have the same number of zeros within $C$.

In our case let the contour $C$ be composed of the unit circle $\Gamma$, (centre at the origin), indented at the point $(1, 0)$ by a small semi-circle $\gamma$ whose radius is $\epsilon$ say.

<figure><img src="image_2.png" /></figure>

Further let us take $f(z) = z$, $g(z) = -\Phi(z)$

On $\Gamma$, $|f(z)| = 1$

and $|g(z)| = |\Phi(z)| = |\sum p(x) z^x| = |\sum p(x) e^{i\theta x}| = |\psi(\theta)| < 1$

i.e on $\Gamma$

$$|g(z)| < 1$$

Again let the radius vector of any point on $\gamma$ by $r$, where

$$0 < r < 1$$

On $\gamma$, $|f(z)| = r$

and $|g(z)| = |\Phi(z)| \le \sum p(x) r^x = \Phi(r)$

We shall show now that on $\gamma$

$$r > \Phi(r)$$

To see this consider the function $\delta(r)$ given by

$$\delta(r) = r - \Phi(r)$$

Now $\delta(1) = 0$

and $\delta'(1) = 1 - \mu < 0$, where $\mu$ is the mean input.
Hence $\delta(r)$ is a decreasing function in the neighbourhood of $r=1$. Thus for values of $r$ slightly less than one, $\delta(r)$ is positive, and since $\epsilon$ is arbitrarily small, we can always have

$$r > \Phi(r)$$

by choosing $\epsilon$ sufficiently small.

Thus we have shown that on all the contour C, the condition $|g(z)| < |f(z)|$ is satisfied.

Moreover the functions $f(z)$ and $g(z)$ are regular within and on C, also $f(z)$ has no zeros on C.

Hence equation (5.4) has only one root inside the unit circle.

This root is real, and to see this consider $\delta(r)$

$$\text{We have} \qquad \delta(1) = 0$$

The function $\delta(r)$ is positive for values of $r$ slightly less than one, but $\delta(0) < 0$

$$\text{Hence} \qquad \delta(r) = 0$$

must have a root between $r=0$, and $r=1$.

Let this root be denoted by $\alpha$. Thus

$$P(t) = 1 - B \alpha^{t+1} \qquad (5.7)$$

To obtain $B$, we observe that this solution must satisfy the equation

$$P(o) = \int_{1}^{\infty} P(y - 1) \, dF(y) \qquad (5.8)$$

This gives $B = 1$, and we thus obtain (5.6).

**VI.** Finally, the merit of this general formulation given by equation (1.2) could be seen by working out an example in the theory of provisioning. This is the case of a "shop" which is receiving regular supplies, but the sale $y$ is a random variable with a distribution $F(y)$.

We take the constant supply as unity.

$$\text{Hence} \qquad u = y - 1$$

$$\text{with the condition} \qquad E(y) < 1$$

$$\text{and} \quad G(u) = F(u + 1)$$

Equation (4.4) takes the form

$$P(t) = \int_{u \le t} P(t - u) \, dF(u + 1)$$
Hence for $t \ge 1$, we have

$$ 
\begin{aligned}
P(t - 1) &= \int_{u+1 \le t} P(t - u - 1) \, dF(u + 1) \\
&= \int_{y \le t} P(t - y) \, dF(y)
\end{aligned}
$$

Since $F(y)$ is zero for negative values of $y$, we finally obtain

$$ P(t - 1) = \int_{0}^{t} P(t - y) \, dF(y), \quad t \ge 1 \quad (6.1) $$

In the case of storing articles which are counted by units, the sale distribution would necessarily be discrete. In this case (6.1) takes the form

$$ P(t - 1) = \sum_{y=0}^{t} P(t - y) f(y) \quad (6.2) $$

where $f(y)$ is the p.d.f. of $y$. Let $\Phi(z)$, $q(z)$ be the probability generating functions of $y$, $t$ i.e

$$ \Phi(z) = \sum_{y=0}^{\infty} z^y f(y) , \quad q(z) = \sum_{t=0}^{\infty} z^t P(t) , \quad |z| < 1 \quad (6.3) $$

Multiplying both sides of (6.2) by $z^t$ and summing from $t=1$ to $\infty$, we get

$$ zq(z) = q(z) \Phi(z) - f(o) P(o) $$

or

$$ q(z) = \frac{f(o) P(o)}{\Phi(z) - z} \quad (6.4) $$

This is a general equation of all discrete distributions. It is obvious that $P(t)$ is the coefficient of $z^t$ in the expansion of R.H.S of (6.4) as a power series of $z$.

For example, if the sale distribution is taken as Poisson's

$$ \text{i.e.} \qquad f(y) = \frac{e^{-\lambda} \lambda^y}{y!} , \quad y = 0, 1, 2, \dots $$

The condition $\lambda < 1$ is necessary for the existence of a stationary distribution,

$$ \text{In this case} \qquad \Phi(z) = e^{(z-1)} $$

$$ \text{Thus} \qquad q(z) = e^{-\lambda} P(o) \sum_{r=0}^{\infty} z^r e^{-\lambda (r+1)} (z-1) $$

and we finally obtain

$$ P(t) = P(o) e^{-\lambda} \sum_{j=0}^{t} \{ \lambda j - \lambda(t+1) \}^{j} \frac{e^{-\{\lambda j - \lambda(t+1)\}}}{j!} \quad (6.5) $$
get

On using Jensen's theorem, and letting $t \to \infty$ in (6.5) we

$$P(o) = (1 - \lambda) e^{\lambda} \qquad (6.6)$$

Hence

$$P(t - 1) = (1 - \lambda) \sum_{j=0}^{t-1} \{\lambda j - \lambda t\}^{j} \frac{e^{-\{\lambda j - \lambda t\}}}{j!}, \quad t > 1 \qquad (6.7)$$

Finally we observe that R.H.S. of (6.6) is less than one as long as $0 < \lambda < 1$.

## REFERENCES

[1] ANIS (A.A.). - The variance of the maximum of partial sums of a finite number of independent normal variates. *Biometrika*, **42**, (1955), 96.

[2] SIMAIKA (J.B.). - On the problem of over-year storage, *Proc. Egyptian Acad. Sci.*, **5**, (1949), 43.

[3] WEESAKUL (B.). - The explicit time-dependent solution for a finite dam with geometric inputs. *Aust. Math. Soc. Summer Research Institute Report (1961)*, (see also *Ann. Math. Statist.*, **32**, (1961), 765-769.

[4] LINDLEY (D.V.). - The theory of queues with a single server, *Proc. camb. Phil. Soc.*, **48**, (1952), 277-289.

[5] SMITHIES (F.). - *Proc. Lond. Math. Soc.*, **46** (1940), 409.

[6] BROCKMEYER (E.), HALSTROM (H.L.) and JENSEN (A.). - *The life and works of A.K. Er lang* (1948).

<figure><img src="image_3.png" /></figure>

Reçu le 3 mars 1965

Ministère du Trésor

Lazoghly, Le Caire

République Arabe Unie

IMPRIMERIE LOUIS-JEAN

Ouvrages scientifiques

TYPO-OFFSET

GAP (Hautes-Alpes)

Dépôt légal n° 132

1966