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 with
rate (λ)(λ). 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 soon
as available, if there is only one. The servers operate independently (of each
other and the arrival process), each with exponential service time. We suppose
each server has the same distribution, exponential (μ)(μ). Such a system may
be analyzed as a Markov birthdeath process. An analysis of the longrun
probabilities 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 and
twoserver systems are utilized to make these calculations quickly and to
present them in a useful way.
Notation
 Nt=Nt= number in system (in service and waiting) at time t
 Qt=Qt= number waiting to be served at time t
 πj=limt→∞pij(t)=πj=limt→∞pij(t)= longrun probability of being in state j
 Wt=Wt= waiting time for service for customer who arrives at
time t
 Dt=Dt= waiting time plus service time for customer who arrives
at time t
 A=A= random variable with distribution of interarrival times
 S=S= random variable with distribution of service times
Longrun probabilities πj=P(Nt=j)πj=P(Nt=j) for large t,st,s servers,
E[A]=1/λ,E[S]=1/μE[A]=1/λ,E[S]=1/μ
For s=1s=1,

ρ
=
E
[
S
]
/
E
[
A
]
=
λ
/
μ
ρ
=
E
[
S
]
/
E
[
A
]
=
λ
/
μ

π
0
=
1

ρ
π
n
=
(
1

ρ
)
ρ
n
π
0
=
1

ρ
π
n
=
(
1

ρ
)
ρ
n

N
t
is
approximately
geometric
(
1

ρ
)
N
t
is
approximately
geometric
(
1

ρ
)
For s>1s>1,

ρ
=
E
[
S
]
/
s
E
[
A
]
=
λ
/
s
μ
ρ
=
E
[
S
]
/
s
E
[
A
]
=
λ
/
s
μ

π
n
=
π
0
(
s
ρ
)
n
/
n
!
=
π
0
(
λ
/
μ
)
n
/
n
!
0
≤
n
≤
s
π
0
(
s
s
/
s
!
)
ρ
n
=
π
0
[
(
s
ρ
)
s
/
s
!
]
ρ
n

s
s
<
n
π
n
=
π
0
(
s
ρ
)
n
/
n
!
=
π
0
(
λ
/
μ
)
n
/
n
!
0
≤
n
≤
s
π
0
(
s
s
/
s
!
)
ρ
n
=
π
0
[
(
s
ρ
)
s
/
s
!
]
ρ
n

s
s
<
n
For s=2s=2

π
0
=
1

ρ
1
+
ρ
=
2
μ

λ
2
μ
+
λ
π
0
=
1

ρ
1
+
ρ
=
2
μ

λ
2
μ
+
λ
For s=3s=3

π
0
=
1

ρ
1
+
2
ρ
+
3
2
ρ
2
π
0
=
1

ρ
1
+
2
ρ
+
3
2
ρ
2
For s=4s=4

π
0
=
1

ρ
1
+
3
ρ
+
8
ρ
2
+
8
3
ρ
3
π
0
=
1

ρ
1
+
3
ρ
+
8
ρ
2
+
8
3
ρ
3
For large t, with the system in equilibrium
E
[
D
t
]
=
E
[
A
]
E
[
N
t
]
and
E
[
W
t
]
=
E
[
A
]
E
[
Q
t
]
≈
E
[
S
]
E
[
N
t
]
E
[
D
t
]
=
E
[
A
]
E
[
N
t
]
and
E
[
W
t
]
=
E
[
A
]
E
[
Q
t
]
≈
E
[
S
]
E
[
N
t
]
(1)
For s=1s=1

E
[
N
t
]
=
ρ
1

ρ
=
λ
μ

λ
E
[
N
t
]
=
ρ
1

ρ
=
λ
μ

λ

E
[
Q
t
]
=
ρ
E
[
N
t
]
P
(
N
t
>
0
)
=
ρ
E
[
Q
t
]
=
ρ
E
[
N
t
]
P
(
N
t
>
0
)
=
ρ

E
[
W
t
]
=
E
[
S
]
E
[
N
t
]
=
λ
/
μ
μ

λ
E
[
W
t
]
=
E
[
S
]
E
[
N
t
]
=
λ
/
μ
μ

λ
 D_{t} is approximately exponential (μλ)(μλ)
For s>1s>1

C
=
P
(
W
t
>
0
)
=
π
0
(
s
ρ
)
s
s
!
(
1

ρ
)
=
E
[
Q
t
]
1

ρ
ρ
=
s
μ
(
1

ρ
)
E
[
W
t
]
C
=
P
(
W
t
>
0
)
=
π
0
(
s
ρ
)
s
s
!
(
1

ρ
)
=
E
[
Q
t
]
1

ρ
ρ
=
s
μ
(
1

ρ
)
E
[
W
t
]

P
(
W
t
>
v
)
=
C
e

(
μ
s

λ
)
v
v
≥
0
P
(
W
t
>
v
)
=
C
e

(
μ
s

λ
)
v
v
≥
0
 P(Dt>v)=eμv1+Cμ1e[μ(s1)λ]vμ(s1)λforλ≠μ(s1)P(Dt>v)=eμv1+Cμ1e[μ(s1)λ]vμ(s1)λforλ≠μ(s1)

P
(
D
t
>
v
)
=
e

μ
v
[
1
+
μ
C
v
]
for
λ
=
μ
(
s

1
)
P
(
D
t
>
v
)
=
e

μ
v
[
1
+
μ
C
v
]
for
λ
=
μ
(
s

1
)

E
[
Q
t
]
=
π
0
(
s
ρ
)
s
s
!
ρ
(
1

ρ
)
2
E
[
Q
t
]
=
π
0
(
s
ρ
)
s
s
!
ρ
(
1

ρ
)
2

E
[
N
t
]
≈
E
[
Q
t
]
+
λ
μ
=
E
[
Q
t
]
+
s
ρ
E
[
N
t
]
≈
E
[
Q
t
]
+
λ
μ
=
E
[
Q
t
]
+
s
ρ
L = input('Enter lambda '); % Type desired value, no extra space
M = input('Enter mu '); % Type desired value, no extra space
a = [' lambda mu'];
b = [L M];
disp(a)
disp(b)
r = L/M; % Rho
EN = 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 values
PD = exp(M*(1  r)*v); % Calculates P(Dt > v)
S = [' v P(D>v)'];
s = [v; PD]';
disp(S)
disp(s)
queue1
Enter lambda 0.1
Enter mu 0.2
lambda mu
0.1000 0.2000
rho EN EQ EW ED
0.5000 1.0000 0.5000 5.0000 10.0000
Enter row matrix of values v [8 16 24]
v P(D>v)
8.0000 0.4493
16.0000 0.2019
24.0000 0.0907
Note that the procedure will not calculate P(D>v)P(D>v) if
λ=μλ=μ.
L = input('Enter lambda '); % Type desired value, no extra space
M = input('Enter mu '); % Type desired value, no extra space
a = [' lambda mu'];
b = [L M];
disp(a)
disp(b)
r = L/(2*M);
EQ = (2*r^3)/(1  r^2);
EN = EQ + 2*r;
EW = EQ/L;
ED = EN/L;
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 ');
t = 2*M*EW*(1  r)/(1  2*r);
PD2 = exp(M*v).*(1 + t.*(1  exp(M*v + L*v))); % Calculates P(D > v) for L not equal M
S = [' v P(D>v)'];
s = [v; PD2]';
disp(S)
disp(s)
queue2
Enter lambda 0.1
Enter mu 0.2
lambda mu
0.1000 0.2000
rho EN EQ EW ED
0.2500 0.5333 0.0333 0.3333 5.3333
Enter row matrix of values v [4 8 16]
v P(D>v)
4.0000 0.4790
8.0000 0.2241
16.0000 0.0473
A queueing system has Poisson arrivals, rate λ and exponential (μ)(μ)
service times.
 In system one, there is one server, with expected service time 1/μ=11/μ=1
minute. Determine
[E[N],E[Q],E[W],E[D],andP(D>v),v=1,3,5,10[E[N],E[Q],E[W],E[D],andP(D>v),v=1,3,5,10
(2)
for expected arrival rates λ=0.6,0.9,0.99λ=0.6,0.9,0.99 customers per minute.
 In system two there are two servers, each with expected service time
1/μ=21/μ=2 minutes. Calculate the same quantities as for system one and compare
the results for the two systems.
queue1
Enter lambda 0.6
Enter mu 1
lambda mu
0.6000 1.0000
rho EN EQ EW ED
0.6000 1.5000 0.9000 1.5000 2.5000
Enter row matrix of values v [1 3 5 10]
v P(D>v)
1.0000 0.6703
3.0000 0.3012
5.0000 0.1353
10.0000 0.0183
Ov = ones(1,length(v));
R = r*Ov; % Row vector with all terms = r
r1 = R;
E11 = B;
v11 = PD;
queue1
Enter lambda 0.9
Enter mu 1
lambda mu
0.9000 1.0000
rho EN EQ EW ED
0.9000 9.0000 8.1000 9.0000 10.0000
Enter row matrix of values v v % Calls for previously entered v
v P(D>v)
1.0000 0.9048
3.0000 0.7408
5.0000 0.6065
10.0000 0.3679
R = r*Ov;
r2 = R;
E12 = B;
v12 = PD;
queue1
Enter lambda 0.99
Enter mu 1
lambda mu
0.9900 1.0000
rho EN EQ EW ED
0.9900 99.0000 98.0100 99.0000 100.0000
Enter row matrix of values v v
v P(D>v)
1.0000 0.9900
3.0000 0.9704
5.0000 0.9512
10.0000 0.9048
R = r*Ov;
r3 = R;
E13 = B;
v13 = PD;
queue2 % Begin calculations for second system
Enter lambda 0.6
Enter mu 0.5
lambda mu
0.6000 0.5000
rho EN EQ EW ED
0.6000 1.8750 0.6750 1.1250 3.1250
Enter row matrix of values v v
v P(D>v)
1.0000 0.7501
3.0000 0.3988
5.0000 0.2019
10.0000 0.0328
E21 = B; % Not necessary to determne r1, r2, r3, since
v21 = PD2; % they are the same as for system one.
queue2
Enter lambda 0.9
Enter mu 0.5
lambda mu
0.9000 0.5000
rho EN EQ EW ED
0.9000 9.4737 7.6737 8.5263 10.5263
Enter row matrix of values v v
v P(D>v)
1.0000 0.9245
3.0000 0.7749
5.0000 0.6410
10.0000 0.3916
E22 = B;
v22 = PD2;
queue2
Enter lambda 0.99
Enter mu 0.5
lambda mu
0.9900 0.5000
rho EN EQ EW ED
0.9900 99.4975 97.5175 98.5025 100.5025
Enter row matrix of values v v
v P(D>v)
1.0000 0.9920
3.0000 0.9743
5.0000 0.9557
10.0000 0.9094
E23 = B;
v23 = PD2;
C = [E11; E21; zeros(E11); E12; E22; zeros(E11); E13; E23]; % Zeros are spacers
disp(A)
rho EN EQ EW ED
disp(C)
0.6000 1.5000 0.9000 1.5000 2.5000
0.6000 1.8750 0.6750 1.1250 3.1250
0 0 0 0 0
0.9000 9.0000 8.1000 9.0000 10.0000
0.9000 9.4737 7.6737 8.5263 10.5263
0 0 0 0 0
0.9900 99.0000 98.0100 99.0000 100.0000
0.9900 99.4975 97.5175 98.5025 100.5025
H = [' rho v P(D1>v) P(D2>v)'];
PDV = [r1 r2 r3; v v v; v11 v12 v13; v21 v22 v23]';
disp(H)
rho v P(D1>v) P(D2>v)
disp(PDV)
1.0000 1.0000 0.6703 0.7501
1.0000 3.0000 0.3012 0.3988
1.0000 5.0000 0.1353 0.2019
1.0000 10.0000 0.0183 0.0328
0.9000 1.0000 0.9048 0.9245
0.9000 3.0000 0.7408 0.7749
0.9000 5.0000 0.6065 0.6410
0.9000 10.0000 0.3679 0.3916
0.9900 1.0000 0.9900 0.9920
0.9900 3.0000 0.9704 0.9743
0.9900 5.0000 0.9512 0.9557
0.9900 10.0000 0.9048 0.9094