Thanks for any suggestions.
S
You could try an iterative method: let pi(n) = c(n)*pi(0), where c(0)
= 1. Then if the arrival rate is 'a' and the service rate per server
is 'b' (with a < m*b), then we have c(n) = c(n-1)*(a/b/n) for n <= m
[where u/v/w means u/(v*w)] and c(n) = c(n-1)*(a/b/m) for n > m. We
then have pi(0) = 1/sum_{n=0..infinity} c(n), etc. The probability
that all servers are busy is pi(0)* sum_{n=m .. inf} c(n), which can
be done because it is a geometric series, etc. Similarly, the expected
queue length Lq is obtainable via a computable series of the form sum_
{n=m..inf} (n-m)*s^(n-m), where s = a/b/m. Then you can obtain L, W,
Wq, etc. in the usual way.
R.G. Vickson
Sandy
Why? Roundoff errors should not be much of a problem, as all you are
doing is adding and multiplying positive numbers. Perhaps overflow or
underflow could be an issue, but using double precision or higher
should get around that.
R.G. Vickson
Overflows and underflows can be avoided by doing the iterations on the
Erlang loss-formula (which keeps the iterates in the interval (0,1))
and then relating the Erlang-B formula to the Erlang-C formula. The
details are spelled out in exercises 4-25 and 4-26 of "Stochastic
Models in Operations Research, Vol. 1" by Heyman and Sobel.
Dan Heyman
Some years ago I wrote a technical paper on this theme; the key
message is that you use the values of a, b and m to determine the most
probable state, and calculate all other probabilities relative to
that. So you then scale the probabilities by the sum of the ratios,
all of which are in the range (0,1) Since the ratios decrease from
the most probable state towards 0 and towards \infty, you can find all
the parameters readily enough.
David K Smith
Is it easy to determine the most probable state before computing the
state probabilities? For m and a/b large, I would estimate the most
probable state by a/b because it's the (or perhaps next to the) most
probable state of the M/M/oo queue.
If s = a/(m*b) < 1, pi(n+1) = s*pi(n) < pi(n) for n >= m implies that
the largest pi(n) is at n <= m; and for n < m we have pi(n) = (a/b/n)
*pi(n-1)>= pi(n-1) if n <= a/b, so it seems that the largest pi(n) is
at n = floor(a/b) or floor(a/b)+1.
R.G. Vickson
Thanks
Sandy
> R.G. Vickson- Hide quoted text -
>
> - Show quoted text -
You don't; that is the whole point of the iterative method! Leave p(k)
unevaluated (where p(k) refers to a specific state probability---such
as the most probable value, or p(0) or whatever you choose). You then
have p(j+1) = r(j)*p(j) for j >= k and p(j-1) = p(j)/r(j-1) for j <=
k, where r(j) = a/b/j for j < m and r(j) = a/b/m = s for j >= m. All p-
values thus are expressed as computed multiples of the unknown p(k),
and you then re-scale to get them to sum to one. In fact, you might as
well set p_old(k) = 1, then get all the other p_old(j) by the
iteration. Let S = sum{p_old(j): j = 0..infinity}. This can be
computed during the computations of the p_old(j) themselves.
Furthermore, S = sum{p_old(j):j=0..m-1} + p_old(m)(1 + s + s^2 + ...)
= sum{p_old(j):j=0..m-1} + p(m)/(1-s). Then set p(j) = p_old(j)/S for
all j. These p(j) sum to 1 and satisfy the steady-state equations. All
the other quantities of interest can be computed in terms of the scale
factor S and the p_old(j). During the iterative computation, if p_old
(j-1) = p_old(j)/r(j-1) leads to p_old(j-1) in the range of underflow,
just stop the computations; to within machine accuracy you can set the
remaining p's to zero. Of course, you don't need numerical values of p
(j) for j >> m because you know that p(j) = p(m)*s^(j-m) for j > m.
Most of this was explained in my first posting, although there I
assumed
you start with an unknown value of p(0) rather than an unknown p(k).
R.G. Vickson