# Distribution of waiting time of "final" customer in finite capacity $M/M/2$ queue with $\mu_1 = 1, \mu_2 = 2, \lambda = 3$

A queue with a maximum capacity of $$5$$ customers has $$2$$ servers operating at rates $$\mu_1 = 1$$ customer/hour and $$\mu_2 = 2$$ customers/hour respectively. Service times are exponentially distributed. Moreover, assume that all customers are served on a first-come, first-served (FCFS) basis. Customers arrive according to a Poisson process at rate $$λ = 3$$ customers/hour. Suppose that there are $$4$$ customers in the salon ($$2$$ being served, $$2$$ waiting) and new customer arrives. What is the distribution of the waiting time of this new customer (the time spent until the customer is served)?

I do not know how to account for the fact that the two servers have different rates. If it were only $$1$$ server with $$3$$ other waiting customers, then the distribution would just be the sum of four independent $$Exp(\mu)$$ variables, which is just an $$Erlang(4, \mu)$$ variable. Also, the question kind of doesn't make sense to me because we are not told how long the $$2$$ customers who are currently being served have been there.

Since the question did not specify, I will assume that customers arriving to the system when there are multiple available servers have equal probability of being served by each available server. We can model this system with a continuous-time Markov chain. Let $$\lambda$$ be the arrival rate of customers to the system, $$\mu_A$$ the service rate of the slower server, and $$\mu_B$$ the service rate of the faster server. Let $$S = \{0, 1_A, 1_B, 2, 3, 4, 5\}$$ and define transition rates $$q_{ij}$$ for each pair of states $$i$$ and $$j$$ by $$\gamma_{ij} = \begin{cases} 0,& i=j\\ \frac\lambda2,& (i,j) \in \{0\}\times\{1_A,1_B\}\\ \lambda,& (i,j)\in \{1_A,1_B\}\times \{2\}\\ \lambda,& i\in\{2,3,4\} \text{ and } j=i+1\\ \mu_A,& (i,j) \in\{(1_A,0), (2,1_B)\}\\ \mu_B,& (i,j) \in\{ (1_B,0), (2,1_A)\}\\ \mu_A+\mu_B,& i\in\{3,4,5\}\text{ and } j=i-1. \end{cases}$$ Define a $$|S|\times|S|$$ matrix $$Q$$ by $$Q_{ij} = \begin{cases} \gamma_{ij},& i\ne j\\ -\sum_{j\in S\setminus\{i\}} \gamma_{ij},& i=j. \end{cases}$$ We call $$Q$$ the transition rate matrix of the process, as we shall see shortly.
Let $$Z_t$$ be the state of the system at time $$t$$, then $$\{Z_t: t\geqslant 0\}$$ is a continuous-time Markov chain; write $$P(t)$$ for the matrix with entries $$p_{ij} = \mathbb P(Z_t = j\mid Z_0=i)$$. Then $$P(t)$$ satisfies the backward Kolmogorov equation $$P'(t) = QP(t)\tag 1$$ with initial condition $$P'(0) = Q$$. From general theory of differential equations we know that there exists a unique solution to $$(1)$$, namely $$P(t) = e^{Qt}$$, where $$e^{Qt} = \sum_{n=0}^\infty \frac{t^nQ^n}{n!}$$ denotes the matrix exponential function. It is a bit unwieldy to write $$P(t)$$ explicitly, but we need not do so to answer the question at hand. Set $$T_0=0$$ and for $$n\geqslant$$ let $$T_n = \inf\{t>T_{n-1} : Z_t\ne Z_{T_{n-1}}\}$$ be the jump times of the process. We can define an embedded Markov chain $$\{X_n:n=0,1,\ldots\}$$ by setting $$X_n = Z_{T_n}$$. Let $$R$$ be the transition matrix of the embedded chain. Recall that if $$U\sim\mathsf{Expo}(\alpha)$$ and $$V\sim\mathsf{Expo}(\beta)$$ are independent, then $$\mathbb P(U (this can be verified by integrating over the joint density of $$(U,V)$$). It follows then that $$R_{ij} = \begin{cases} 0,& i=j\\ \frac{\gamma_{ij}}{\sum_{k\in S}\gamma_{ik}},& i\ne j. \end{cases}$$ In our example, $$R = \left( \begin{array}{ccccccc} 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 \\ \frac{\mu _A}{\mu _A+\lambda } & 0 & \frac{\lambda }{\mu _A+\lambda } & 0 & 0 & 0 & 0 \\ \frac{\mu _B}{\mu _B+\lambda } & 0 & 0 & \frac{\lambda }{\mu _B+\lambda } & 0 & 0 & 0 \\ 0 & \frac{\mu _B}{\mu _A+\mu _B+\lambda } & \frac{\mu _A}{\mu _A+\mu _B+\lambda } & 0 & \frac{\lambda }{\mu _A+\mu _B+\lambda } & 0 & 0 \\ 0 & 0 & 0 & \frac{\mu _A+\mu _B}{\mu _A+\mu _B+\lambda } & 0 & \frac{\lambda }{\lambda+\mu _A+\mu _B} & 0 \\ 0 & 0 & 0 & 0 & \frac{\mu _A+\mu _B}{\mu _A+\mu _B+\lambda } & 0 & \frac{\lambda }{\mu _A+\mu _B+\lambda } \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ \end{array} \right)$$ We wish to find the distribution of $$\tau = \inf\{t>0: Z_t \in\{1_A,1_B\}\mid Z_0 = 5\}.$$
This is not an easy task. Consider the simpler question of finding the distribution of $$N= \inf\{n>0: X_n\in\{1_A,1_B\}\mid X_0=5\}.$$ This is the number of jumps (services and arrivals) until one of the servers becomes available to serve the new customer. $$N$$ has a discrete phase-type distribution, i.e. it is the distribution of the time hitting time of the absorbing state in a terminating Markov chain (one in which all but the absorbing state are transient). For this terminating Markov chain, lump together the states $$1_A$$ and $$1_B$$ into the state $$1$$, and reorder the states as $$(5,4,3,2,1)$$. Now consider the substochastic matrix of $$R$$, call it $$R_a$$, obtained by this lumping and reordering process (as well as removing the state $$0$$, as it is not relevant). Then $$R_a = \left( \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \\ \frac{\lambda }{\mu _A+\mu _B+\lambda } & 0 & \frac{\mu _A+\mu _B}{\mu _A+\mu _B+\lambda } & 0 & 0 \\ 0 & \frac{\lambda }{\mu _A+\mu _B+\lambda } & 0 & \frac{\mu _A+\mu _B}{\mu _A+\mu _B+\lambda } & 0 \\ 0 & 0 & \frac{\lambda }{\mu _A+\mu _B+\lambda } & 0 & \frac{\mu _A+\mu _B}{\mu _A+\mu _B+\lambda } \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right).$$ We can write $$R_a = \begin{pmatrix} T & \mathbf T^0\\ \mathbf 0 & \mathbf 1, \end{pmatrix}$$ where $$T$$ is the substochastic matrix corresponding to transitions between transient states and $$\mathbf T^0 + T\mathbf 1=\mathbf 1$$. It can be shown that $$p_n:=\mathbb P(N=n)$$ is given by $$p_n = \tau T^{n-1}\mathbf T^0,$$ where $$\tau$$ is the initial distribution (in our case, $$\tau = \begin{pmatrix}1&0&0&0&0\end{pmatrix}^T$$). There does not appear to be a concise closed form for the $$n^{\mathrm{th}}$$ power of $$T$$, although we can compute $$\mathbb E[N]$$ without much trouble. Define $$N_a := \sum_{k=0}^\infty T^k = (\mathbf{I_4}-T)^{-1},$$ where $$\mathbf{I_4}$$ is the $$4\times 4$$ identity matrix. (Note that this series converges because there is a row which sums to strictly less than one.) Then the expected number of steps until absorption when starting in state $$5$$ is the first entry of the vector $$N\mathbf 1$$. Doing the computations, we find that $$\mathbb E[N] = \frac{2 \left(\mu _A+\mu _B+\lambda \right) \left(\left(\mu _A+\mu _B\right) \left(2 \mu _A+2 \mu _B+\lambda \right)+\lambda ^2\right)}{\left(\mu _A+\mu _B\right){}^3},$$ which is equal to $$16$$ when $$\lambda=3$$, $$\mu_A=1$$, and $$\mu_B=2$$.
Unfortunately, even the exact distribution of $$N$$ would not lend much insight as to the distribution of $$\tau$$. This is because the holding time in state $$5$$ is exponential with rate $$\mu_A+\mu_B$$, but the holding time in states $$4$$ and $$3$$ have a hyperexponential distribution with density $$f(t) = \frac{\lambda(\mu_A+\mu_B)}{\lambda+\mu_A+\mu_B}\left(e^{-(\mu_A+\mu_B)t} + e^{-\lambda t} \right)$$ and mean $$\frac{2 \mu _A \mu _B+\mu _A^2+\mu _B^2+\lambda ^2}{\lambda \left(\mu _A+\mu _B\right) \left(\mu _A+\mu _B+\lambda \right)}.$$ Although your question concerns the transient behavior of the process $$\{Z_t\}$$, it may be worthwhile to compute the stationary distribution (the rows of $$\lim_{t\to\infty} e^{Qt}$$), and proceed from there. To my knowledge there is no general way to find the distribution of the hitting time of a state in a continuous-time Markov chain (but if there is, then feel free to correct me!)
• Thank you for the detailed answer. I managed to solve it in a different way. If it were one server, then the question was easy. So I 'combined' the two servers into one working at rate $= 1 + 2 = 3$ and then used the fact that this 'combined' server had to serve exactly $3$ customers before serving the last one. Thus the distribution of the waiting time is $Erlang(3, 3)$.