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.
 A: I was not able to find an analytical solution to this question, but I hope the following analysis is useful in your attempts.
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<V) = \frac{\alpha}{\alpha+\beta}$ (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!)
