Distribution of server utilisations in an M/M/c queuing model with an unusual dispatching discipline I'm studying an M/M/c queuing model with an unusual (?) dispatching discipline:


*

*Servers are numbered 1...c

*The servers have an identical mean service time, exponentially distributed (as usual), which does not vary with time or load

*If all servers are busy, the transaction is allocated to the first server that becomes free

*If any servers are free, the transaction is allocated to the free server with the lowest number.


I am especially interested in the mean server utilisation ($\rho_i$) for each server (which could be derived from the proportion $p_i$ of jobs served by server $i$), and also its distribution (though I guess that is more difficult).
What results are available which give the distribution of traffic going to each server? 

This is (an attempt at) a partial answer.
I think that you should be able to decouple the server selection aspect from the composite load aspect, since the servers are identical with exponentially distributed service time $\mu$.  That is, one can analyze the number in system (irrespective of distribution across the servers) as an ordinary M/M/$c$ system with the usual distribution:
$$
p_k = \begin{cases}
\hfill p_0 \frac{\sigma^k}{k!} \hfill & k \leq c \\
\hfill p_0 \frac{\sigma^k}{c!c^{k-c}} \hfill & k \geq c
\end{cases}
$$
where $\sigma = \lambda/\mu$ and
$$
p_0 = \left(\sum_{k=0}^{c-1} \frac{\sigma^k}{k!}
          + \sum_{k=c}^\infty \frac{\sigma^k}{c!c^{k-c}}\right)^{-1}
    = \left(\frac{c\sigma^c}{c!(c-\sigma)}
          + \sum_{k=0}^{c-1} \frac{\sigma^k}{k!}\right)^{-1}
$$
For the simplest case $c=2$, we can then disentangle the individual states $(1, 0)$ (server $1$ busy, server $2$ idle) and $(0, 1)$ (server $1$ idle, server $2$ busy) by writing
$$
\lambda p_0 = \mu p_{1, 0} + \mu p_{0, 1} = \mu p_1
$$
$$
(\lambda+\mu) p_{1, 0} = \lambda p_0 + \mu p_2
$$
$$
(\lambda+\mu) p_{0, 1} = \mu p_2
$$
$$
2\mu p_2 = \lambda p_{1, 0} + \lambda p_{0, 1} = \lambda p_1
$$
Together, these equations yield
$$
p_{1, 0} = \frac{2+\sigma}{2+2\sigma} \, p_1
$$
$$
p_{0, 1} = \frac{\sigma}{2+2\sigma} \, p_1
$$
I'm still thinking about (a) whether this is all correct for $c = 2$, and (b) if it is, how one might generalize for all $c$.
