G/M/1 Queue with Arrivals Following an Erlang-2 Distribution I look at the following queueing system:
"We have two identical M/M/1 queueing systems, S1 and S2. Jobs arrive according to a Poisson process with rate 2 per hour, and service times are exponential with mean $1/4$ hour. The arrival stream is split into two (stchastically equal) arrival streams, which form the arrival processes for S1 and S2.
Suppose that the spliting is done such that the jobs are alternatingly routed to S1 and S2; as a result te interarrival times at S1 have an Erlang-2 distribution."
Firstly, I want to compute $\rho_1$, the fraction of time the server (at S1) is busy. I know that $\rho_1=\lambda_1 \cdot \mathbb{E}[B]$, where I denote with $\lambda_1$ the arrival rate at S1.
I also know that the interarrival time of the whole system has mean $\mathbb{E}[A]=\frac{1}{\lambda}$, which would give $\mathbb{E}[A]=\frac{1}{2}$. So my idea was, knowing that the mean of the Erlang-2 distribution is $\frac{2}{\lambda_1}$ to compute $\lambda_1$ with
$$\mathbb{E}[A]=\frac{1}{2}=\frac{2}{\lambda_1}+\frac{2}{\lambda_2},$$
where $\lambda_2$ belongs to S2. Now since the situation in S2 is symmetric to the one in S1, I have $\lambda_1=\lambda_2$. That would give me $\lambda=8$ and consequently an unstable system because $\rho=4 \cdot \frac{1}{4}$.
Can somebody help me understand the system? How do I derive the correct arrival rate?
 A: We can consider each station as a $G/M/1$ queue with arrival distribution $\mathrm{Erlang}(2,\lambda)$, that is, the interarrival times have density
$$
f_A(t)=\lambda(\lambda t)e^{-\lambda t}\mathsf 1_{(0,\infty)}(t).
$$
Let $X_0=0$ and $X_n$ be the number of customers in the system just before the $n^{\mathrm{th}}$ arrival. Then $\{X_n:n=0,1,2,\ldots\}$ is the embedded Markov chain on the nonnegative integers with transition probabilities
$$
\mathbb P(X_{n+1}=j\mid X_n=i) = \begin{cases}
p_{i,0},& j=0\\
\beta_{i-j+1},& j>0.
\end{cases}
$$
Here $\beta_j$ is the probability of serving $j$ customers during an interarrival time given that the server remains busy during this interval. Let $T$ be the interarrival time, then conditioned on $\{T=t\}$ the number of customers served $N$ has Poisson distribution with parameter $\mu t$. Thus
\begin{align}
\beta_j &= \int_0^\infty \mathbb P(N=j\mid T=t)f_A(t)\ \mathsf dt\\
&= \int_0^\infty \frac{(\mu t)^j}{j!} e^{-\mu t}\lambda(\lambda t)e^{-\lambda t}\ \mathsf dt\\
&= (j+1)\left(\frac\lambda\mu\right)^2\left(\frac\mu{\lambda+\mu}\right)^{j+2},\ j=0,1,2,\ldots.
\end{align}
Since the rows of the transition matrix should sum to one, it follows that
\begin{align}
p_{i,0} &= 1 - \sum_{j=0}^i \beta_j\\
&= \sum_{i=j+1}^\infty \beta_j\\
&= \sum_{i=j+1}^\infty (i+1)\left(\frac\lambda\mu\right)^2\left(\frac\mu{\lambda+\mu}\right)^{i+2}\\
&= \frac{\mu(\mu+\lambda(i+2))}{(\lambda+\mu)^2}\left(\frac\mu{\lambda+\mu}\right)^i,\ i=0,1,2,\ldots
\end{align}
The stationary distribution $\pi$ satisfies the balance equations
\begin{align}
\pi_0 &= \sum_{i=0}^\infty p_{i,0}\pi_i\\
\pi_n &= \sum_{i=0}^\infty \beta_i\pi_{n-i+1},\ n=1,2,\ldots\tag1.
\end{align}
It is known in the literature that $\pi$ has a geometric distribution, i.e. $\pi_n = (1-\sigma)\sigma^n$ for some $\sigma\in(0,1)$. Substituting this into $(1)$ and dividing by the common factor $\sigma^{n-1}$ yields
\begin{align}
\sigma &= \sum_{i=0}^\infty \sigma^i\beta_i\\
&= \sum_{i=0}^\infty \sigma^i (i+1)\left(\frac\lambda\mu\right)^2\left(\frac\mu{\lambda+\mu}\right)^{i+2}\\
&= \frac{\lambda ^2}{(\lambda + \mu(1-\sigma))^2}.
\end{align}
Solving for $\sigma$ yields the roots
\begin{align}
\sigma &= \frac{\mu ^{3/2} \left(-\sqrt{4 \lambda +\mu }\right)+2 \lambda  \mu +\mu ^2}{2 \mu ^2}\tag2\\\
\sigma &= \frac{\mu ^{3/2} \sqrt{4 \lambda +\mu }+2 \lambda  \mu +\mu ^2}{2 \mu ^2}\tag3.
\end{align}
Since $\sum_{i=0}^\infty\pi_i=1$, it turns out that $(2)$ is the correct root. So we have
$$
\pi_n = \left(\frac{\sqrt{\mu } \sqrt{4 \lambda +\mu }-2 \lambda +\mu }{2 \mu }\right) \left(\frac{2 \lambda +\mu -\sqrt{\mu(1+4\lambda)}  }{2\mu }\right)^n,\ n=0,1,2,\ldots
$$
The fraction of time that the server is busy is given by
$$
1 - \pi_0 = \frac{2 \lambda +\mu -\sqrt{\mu(1+4\lambda)}  }{2\mu }.
$$
