M/M/2/2 interdeparture distribution

Here’s a short derivation of the M/M/2/2 interdeparture distribution that I recently did as a homework exercise. It was a bit hard to find this result online, but it’s a pretty short derivation so I thought it would be helpful to keep around.

In [3]:

import sympy

Derivation of interdeparture distribution for $M/M/2/2$

In [4]:

lambd = sympy.Symbol("lambda")
mu = sympy.Symbol("mu")
s = sympy.Symbol("s")

To derive the departure distribution, we think about the following. Suppose a departure has occurred.

Then either the transition $2 \rightarrow 1$ happened or the transition $1 \rightarrow 0$ happened.

Denote the probability of those events by $\tilde p_1$ and $\tilde p_0$, respectively.

Let $\phi_0(s)$ and $\phi_1(s)$ be the Laplace transforms (MGF) of the time till the next departure starting from states 0 and states 1, respectively.

If we are in state $0$, we have to wait for the next arrival and then we are in state $1$.
If we are in state $1$, two things could happen:
- A departure occurs (the desired event).
- An arrival occurs sending us to $2$, so then we would wait for the next departure from here.

Then, the Laplace transform of the interdeparture distribution is the mixture of these two: $$ \phi(s) = \tilde p_0 \phi_0(s) + \tilde p_1 \phi_1(s) $$

For $\phi_1$:

We have to wait for event (either arrival or departure to occur), which is exponentially distributed with rate $\lambda + \mu$
The probability that event is an arrival is $\lambda/(\lambda + \mu)$. That sends us to state 2, where departure time is Exp with rate $2\mu$
The probability that event is a departure is $\mu/(\lambda + \mu)$. Then the expected waiting time is 0, which for the Laplace transform is $e^0 = 1$.

Sum of R.Vs is multiplication of Laplace transforms. We'll use $$ \frac{\lambda}{\lambda + s} $$ as the Laplace transform of the exponential PDF with rate $\lambda$ (sympy likes the +s version).

In [7]:

phi_1 = (lambd + mu)/(lambd + mu + s)
phi_1 *= ( (mu/(mu+lambd)) + (lambd/(mu+lambd))*(2*mu / (2*mu+s)))
phi_1

Out[7]:

$\displaystyle \frac{\left(\lambda + \mu\right) \left(\frac{2 \lambda \mu}{\left(\lambda + \mu\right) \left(2 \mu + s\right)} + \frac{\mu}{\lambda + \mu}\right)}{\lambda + \mu + s}$

For $\phi_0$:

We wait for the arrival to occur, which is exponentially distributed with rate $\lambda$, then we are state 1 so the remaining time is given by $\phi_1$

In [8]:

phi_0 = (lambd/(lambd+s))*phi_1
phi_0

Out[8]:

$\displaystyle \frac{\lambda \left(\lambda + \mu\right) \left(\frac{2 \lambda \mu}{\left(\lambda + \mu\right) \left(2 \mu + s\right)} + \frac{\mu}{\lambda + \mu}\right)}{\left(\lambda + s\right) \left(\lambda + \mu + s\right)}$

To derive the $\tilde p_0$ and $\tilde p_1$, in steady-state the probability of observing a departure and landing in a certain state is given by

$$ \frac{\text{Departure rate into i}}{\text{Total departure rate}} $$

In steady state: Total departure rate = Total arrival rate, so the denominator is $$ \lambda (1 - \text{Blocking probability}) $$

The blocking is $p_2$, the steady-state time average fraction the system is in $2$, which we can calculate using the Erlang loss formula (derived via CTMC analysis). It's also called Erlang's B formula (https://en.wikipedia.org/wiki/Erlang_(unit)#Erlang_B_formula)

To get the departure rate into $1$: $$ (\text{Fraction of time in 2})(\text{Departure rate from } 2 \rightarrow 1) = p_2 (2\mu) $$

And similarly for into $0$: $$ (\text{Fraction of time in 1})(\text{Departure rate from } 1 \rightarrow 0) = p_1 \mu $$

In [9]:

rho = lambd/mu
pi_0 = (1 + rho + rho**2 / 2)**-1
pi_1 = rho*pi_0
pi_2 = rho * pi_1 / 2

In [10]:

p0 = (pi_1 * mu)/((1-pi_2)*lambd)
p1 = (pi_2 * 2 * mu) / ((1 - pi_2)*lambd)

In [11]:

phi = p0*phi_0 + p1*phi_1
phi = sympy.simplify(phi)

In [12]:

phi

Out[12]:

$\displaystyle \frac{\lambda \mu \left(2 \lambda + 2 \mu + s\right)}{\left(\lambda + \mu\right) \left(\lambda + s\right) \left(2 \mu + s\right)}$

Now we just need to invert the Laplace transform, which sympy has a function for:

In [13]:

from sympy import inverse_laplace_transform
from sympy.abc import t

In [14]:

f = inverse_laplace_transform(phi, s, t)

In [15]:

f = sympy.simplify(f)

In [16]:

Out[16]:

$\displaystyle \frac{\lambda \mu \left(2 \lambda e^{\lambda t} - \left(\lambda + 2 \mu\right) e^{2 \mu t}\right) e^{- t \left(\lambda + 2 \mu\right)} \theta\left(t\right)}{\left(\lambda - 2 \mu\right) \left(\lambda + \mu\right)}$

$\theta(t)$ is the Heaviside function which plays the role of the indicator function, saying that this function is only positive on the positive reals and zero elsewhere.

See: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.delta_functions.Heaviside

In [18]:

sympy.plot(
    f.subs(mu, 1/5).subs(lambd, 1/5),
    (t,0,30)
)

No description has been provided for this image

Out[18]:

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7f53c4958a50>

In [23]:

sympy.srepr(f)

Out[23]:

"Mul(Symbol('lambda'), Symbol('mu'), Pow(Add(Symbol('lambda'), Mul(Integer(-1), Integer(2), Symbol('mu'))), Integer(-1)), Pow(Add(Symbol('lambda'), Symbol('mu')), Integer(-1)), Add(Mul(Integer(2), Symbol('lambda'), exp(Mul(Symbol('lambda'), Symbol('t')))), Mul(Integer(-1), Add(Symbol('lambda'), Mul(Integer(2), Symbol('mu'))), exp(Mul(Integer(2), Symbol('mu'), Symbol('t'))))), exp(Mul(Integer(-1), Symbol('t'), Add(Symbol('lambda'), Mul(Integer(2), Symbol('mu'))))), Heaviside(Symbol('t')))"

We can even get the CDF:

In [20]:

F = sympy.simplify(sympy.integrate(f,t))
F

Out[20]:

$\displaystyle \begin{cases} 0 & \text{for}\: \lambda = 0 \wedge \left(\lambda = 0 \vee \mu = 0\right) \\\frac{\lambda \mu \left(\frac{\lambda}{\mu} - \frac{\lambda e^{- 2 \mu t}}{\mu} - 1 + e^{- \lambda t} - \frac{2 \mu}{\lambda} + \frac{2 \mu e^{- \lambda t}}{\lambda}\right) \theta\left(t\right)}{\left(\lambda - 2 \mu\right) \left(\lambda + \mu\right)} & \text{otherwise} \end{cases}$

In [21]:

sympy.plot(
    F.subs(mu, 1/5).subs(lambd, 1/5),
    (t,0,30)
)

Out[21]:

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7f53c4927110>

In [22]:

sympy.srepr(F)

Out[22]:

"Piecewise(ExprCondPair(Integer(0), And(Equality(Symbol('lambda'), Integer(0)), Or(Equality(Symbol('lambda'), Integer(0)), Equality(Symbol('mu'), Integer(0))))), ExprCondPair(Mul(Symbol('lambda'), Symbol('mu'), Pow(Add(Symbol('lambda'), Mul(Integer(-1), Integer(2), Symbol('mu'))), Integer(-1)), Pow(Add(Symbol('lambda'), Symbol('mu')), Integer(-1)), Add(Mul(Symbol('lambda'), Pow(Symbol('mu'), Integer(-1))), Mul(Integer(-1), Symbol('lambda'), Pow(Symbol('mu'), Integer(-1)), exp(Mul(Integer(-1), Integer(2), Symbol('mu'), Symbol('t')))), Integer(-1), exp(Mul(Integer(-1), Symbol('lambda'), Symbol('t'))), Mul(Integer(-1), Integer(2), Pow(Symbol('lambda'), Integer(-1)), Symbol('mu')), Mul(Integer(2), Pow(Symbol('lambda'), Integer(-1)), Symbol('mu'), exp(Mul(Integer(-1), Symbol('lambda'), Symbol('t'))))), Heaviside(Symbol('t'))), true))"