Basic theory of one and twoserver ques with Poisson arrivals and exponential servers. Matlab calculations provide numerical examples.
A standard model of a queueing system with a single waiting line and one or more
servers assumes that “customers” arrive according to a Poisson process withrate
$\left(\lambda \right)$ . The customer at the head of the line goes to the first
available server, if there are more than one, or to the single server as soonas available, if there is only one. The servers operate independently (of each
other and the arrival process), each with exponential service time. We supposeeach server has the same distribution, exponential
$\left(\mu \right)$ . Such a system may
be analyzed as a Markov birthdeath process. An analysis of the longrunprobabilities and expectations of various quantities after the system has
settled down to equilibrium yields the results below.
Calculation of these quantities is straightforward, but somewhat tedious if
various cases are considered. Matlab procedures for singleserver andtwoserver systems are utilized to make these calculations quickly and to
present them in a useful way.
Notation

${N}_{t}=$ number in system (in service and waiting) at time
t

${Q}_{t}=$ number waiting to be served at time
t

${\pi}_{j}=\underset{t\to \infty}{lim}{p}_{i}j\left(t\right)=$ longrun probability of being in state
j

${W}_{t}=$ waiting time for service for customer who arrives at
time
t

${D}_{t}=$ waiting time plus service time for customer who arrives
at time
t

$A=$ random variable with distribution of interarrival times

$S=$ random variable with distribution of service times
Longrun probabilities
${\pi}_{j}=P({N}_{t}=j)$ for large
$t,s$ servers,
$E\left[A\right]=1/\lambda ,E\left[S\right]=1/\mu $
For
$s=1$ ,

$\rho =E\left[S\right]/E\left[A\right]=\lambda /\mu $

${\pi}_{0}=1\rho {\pi}_{n}=(1\rho ){\rho}^{n}$

${N}_{t}\phantom{\rule{4pt}{0ex}}\text{is}\phantom{\rule{4pt}{0ex}}\text{approximately}\phantom{\rule{4pt}{0ex}}\text{geometric}(1\rho )$
For
$s>1$ ,

$\rho =E\left[S\right]/sE\left[A\right]=\lambda /s\mu $

${\pi}_{n}=\left\{\begin{array}{cc}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}{\pi}_{0}\phantom{\rule{0.166667em}{0ex}}{\left(s\rho \right)}^{n}/n!={\pi}_{0}\phantom{\rule{0.166667em}{0ex}}{(\lambda /\mu )}^{n}/n!\hfill & 0\le n\le s\hfill \\ \phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}{\pi}_{0}\phantom{\rule{0.166667em}{0ex}}({s}^{s}/s!){\rho}^{n}={\pi}_{0}\phantom{\rule{0.166667em}{0ex}}[{\left(s\rho \right)}^{s}/s!]{\rho}^{ns}\hfill & s<n\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}\hfill \end{array}\right)$
For
$s=2$

$\pi}_{0}=\frac{1\rho}{1+\rho}=\frac{2\mu \lambda}{2\mu +\lambda$
For
$s=3$

$\pi}_{0}=\frac{1\rho}{1+2\rho +\frac{3}{2}\phantom{\rule{0.166667em}{0ex}}{\rho}^{2}$
For
$s=4$

$\pi}_{0}=\frac{1\rho}{1+3\rho +8{\rho}^{2}+\frac{8}{3}\phantom{\rule{0.166667em}{0ex}}{\rho}^{3}$
For large
t , with the system in equilibrium
$$E\left[{D}_{t}\right]=E\left[A\right]E\left[{N}_{t}\right]\phantom{\rule{4pt}{0ex}}\text{and}\phantom{\rule{4pt}{0ex}}E\left[{W}_{t}\right]=E\left[A\right]E\left[{Q}_{t}\right]\approx E\left[S\right]E\left[{N}_{t}\right]$$
For
$s=1$

$E\left[{N}_{t}\right]=\frac{\rho}{1\rho}=\frac{\lambda}{\mu \lambda}$

$E\left[{Q}_{t}\right]=\rho E\left[{N}_{t}\right]P({N}_{t}>0)=\rho $

$E\left[{W}_{t}\right]=E\left[S\right]E\left[{N}_{t}\right]=\frac{\lambda /\mu}{\mu \lambda}$

D
_{t} is approximately exponential
$(\mu \lambda )$
For
$s>1$

$C=P({W}_{t}>0)={\pi}_{0}\frac{{\left(s\rho \right)}^{s}}{s!(1\rho )}=E\left[{Q}_{t}\right]\frac{1\rho}{\rho}=s\mu (1\rho )E\left[{W}_{t}\right]$

$P({W}_{t}>v)=C{e}^{(\mu s\lambda )v}v\ge 0$

$P({D}_{t}>v)={e}^{\mu v}\left[1,+,C,\mu ,\frac{1{e}^{\left[\mu \right(s1)\lambda ]v}}{\mu (s1)\lambda}\right]\phantom{\rule{4pt}{0ex}}\text{for}\phantom{\rule{4pt}{0ex}}\lambda \ne \mu (s1)$

$P({D}_{t}>v)={e}^{\mu v}[1+\mu Cv]\phantom{\rule{4pt}{0ex}}\text{for}\phantom{\rule{4pt}{0ex}}\lambda =\mu (s1)$

$E\left[{Q}_{t}\right]={\pi}_{0}\frac{{\left(s\rho \right)}^{s}}{s!}\frac{\rho}{{(1\rho )}^{2}}$

$E\left[{N}_{t}\right]\approx E\left[{Q}_{t}\right]+\frac{\lambda}{\mu}=E\left[{Q}_{t}\right]+s\rho$
Matlab calculations for single server queue (in file queue1.m)
L = input('Enter lambda '); % Type desired value, no extra space
M = input('Enter mu '); % Type desired value, no extra spacea = [' lambda mu'];b = [L M];disp(a)
disp(b)r = L/M; % RhoEN = r/(1  r); % E[N]EQ = r*EN; % E[Q]EW = EQ/L; % E[W]ED = EN/L; % E[D]A = [' rho EN EQ EW ED']; % Identifies entries in B
B = [r EN EQ EW ED];
disp(A)disp(B)v = input('Enter row matrix of values v '); % Type matrix of desired valuesPD = exp(M*(1  r)*v); % Calculates P(Dt>v)S = [' v P(D>v)'];s = [v; PD]';disp(S)
disp(s)
queue1Enter lambda 0.1
Enter mu 0.2lambda mu
0.1000 0.2000rho EN EQ EW ED0.5000 1.0000 0.5000 5.0000 10.0000Enter row matrix of values v [8 16 24]v P(D>v)
8.0000 0.449316.0000 0.2019
24.0000 0.0907