(c) Tobias Hossfeld (Aug 2021)
This script and the figures are part of the following book. The book is to be cited whenever the script is used (copyright CC BY-SA 4.0):
Tran-Gia, P. & Hossfeld, T. (2021). Performance Modeling and Analysis of Communication Networks - A Lecture Note. Würzburg University Press. https://doi.org/10.25972/WUP-978-3-95826-153-2
With the time-discrete analysis of GI/GI/1 delay systems, the distribution function of the waiting time can be determined numerically. Any continuous distributions can be approximated by discrete distributions. Using the example of the M/GI/1 waiting system, the exact waiting time (using Laplace transformation) is compared with the results of the time-discrete analysis. The finer the discretization $ \Delta t $, the more precise the time-discrete analysis. This parameter $ \Delta t $ should be played around in the notebook to see how this affects the accuracy of the numerical results.
The arrival process is a Poisson process with rate $ \lambda = 1 / E[A] $. The service time follows an Erlang-k distribution, $ B \sim E_k(\mu_i) $ with $ E[B] = 1 / \mu $. Each of the $ k $ phases of the Erlang-k distribution follows an exponential distribution with rate $ \mu_i = k \mu = k / E[B] $.
The interarrival time $ A $ and the service time $ B $ are specified. In addition to the parameters of the distributions, the Laplace transform of the service time is also given and implemented as LT_B (s).
from matplotlib import pyplot as plt
import numpy as np
from mpmath import *
import scipy.stats
from discreteTimeAnalysis import *
##################################
# Interarrival time A ~ Exp(lam)
lam = 0.7 # arrival rate
EA = 1/lam # mean interarrival time
cA = 1 # coefficient of variation of exponential distribution
# PDF of A
def at(t, lam=lam):
return lam*np.exp(-lam*t)
# CDF of A
def At(t, lam=lam):
return 1-np.exp(-lam*t)
##################################
# service time B ~ Erlang_k(lam)
EB = 1 # mean service time
mu = 1/EB # service rate
k = 3 # k exponential phases of Erlang-k distribution
cB = 1/np.sqrt(k) # coefficient of variation of the Erlang-k distribution
# Laplace transform of service time B ~ Erlang_k
def LT_B(s, EB=EB, n=k):
return 1/(1+s*EB/n)**n # Erlang-k
# CDF of the service time
def Bt(t, EB=EB, n=k):
rvB = scipy.stats.erlang(n,scale=EB/n) # we use the scipy implementation
return rvB.cdf(t)
##################################
# Stability condition: rho<1
rho =lam/mu
# title of plots
title = f'M/M/1 with $\lambda={lam}, E[B]={EB}$' if k==1 else f'M/E$_{k}$/1 with $\lambda={lam}, E[B]={EB}, cB={cB:.2f}$'
print(title)
M/E$_3$/1 with $\lambda=0.7, E[B]=1, cB=0.58$
For an M/GI/1 system, the Laplace transform of the waiting time is as follows. The Laplace transform of the service time is required for the calculation:
$ LT \{b (t) \} = \Phi_B (s) $.
$ \displaystyle LT\{w(t)\} = \Phi_W(s) = \frac{(1-\rho)\cdot s}{s-\lambda+\lambda \Phi_B(s)}$
# Laplace transform of waiting time
def LT_W(s, lam=lam, EB=EB):
rho = EB*lam
return (1-rho)*s/(s-lam+lam*LT_B(s, EB))
s = np.linspace(0.00001,1,100)
plt.plot(s, LT_B(s), label='$\Phi_B(s)$=LT{b(t)}')
plt.plot(s, LT_W(s), label='$\Phi_W(s)$=LT{w(t)}')
plt.xlabel('s')
plt.ylabel('LaplacetTransform')
plt.legend()
plt.grid()
plt.title(title);
The inverse numerical Laplace transform delivers the probability density function (PDF) of the waiting time
$LT\{w(t)\} = \Phi_W(s)$.
$w(t) = LT^{-1}\{\Phi_W(s)\}$
To calculate the cumulative distribution function (CDF), the inverse Laplace transform of $ \Phi_W (s) / s $ is calculated.
$ \displaystyle W(t) = LT^{-1}\{\frac{\Phi_W(s)}{s}\} $
The mpmath package can be used in Python for the inverse Laplace transform purpose. The function is called invertlaplace.
# PDF of the waiting time : w(t) (inverse numerical Laplace transform)
w = lambda t: invertlaplace(LT_W, t ,method='talbot')
# CDF of the waiting time: LT{W(t)} = LT{w(t)}/s
W = lambda t: invertlaplace(lambda s: LT_W(s)/s, t, method='talbot')
# Comparison with an M/M/1 System: PDF of the waiting time
# for k=1 the systems the M/E_k/1 and M/M/1 system are identical.
def wmm1(t, lam=lam, EB=EB):
mu = 1/EB
rho = lam/mu
return mu*(1-rho)*rho*np.exp(-mu*(1-rho)*t)
# Let's plot it
t_lt = np.linspace(1e-9,10,100)
wt = np.array([w(ti) for ti in t_lt])
plt.plot(t_lt, wt, label='numerical inverse Laplace')
plt.plot(t_lt, wmm1(t_lt), label='M/M/1')
plt.xlabel('waiting time (s)')
plt.ylabel('PDF')
plt.legend()
plt.grid()
plt.title(title);
The probability $ P (A \leq t) = A (t) $ is calculated at the discrete places $ i \Delta t $ for $ i = 0,1, ..., a_\max $. Then the discrete distribution $ A_d $ is specified by the probabilities
$ a_d (i) = P (A_d = i) = A (i \cdot \Delta t) -A ((i-1) \cdot \Delta t) \quad $ for $\; 0 <i <a_\max $.
It is:
$a_d(0)=0 \\ a_d(a_\max) = 1-A\left((a_\max-1)\cdot \Delta t\right).$
#%% discretization constant
deltaT = 0.5
# compute discretized interarrival times and CDF values
amax = 50/deltaT # define maximum value a_max
tdisc = np.arange(0, amax+1)*deltaT # compute A(t) at the point tdisc
A = At(tdisc)
A[-1] = 1 # insert 1 to compute a(a_max)
# compute probabilities (PMF values)
a = np.insert(np.diff(A), 0, 0) # insert 0 to compute a(0)
xa = np.arange(0, amax+1).astype(int)
EAdisc = np.matmul(a,xa) # mean value
# Generate the corresponding distribution object
Adisc = DiscreteDistribution(xa, a)
print(f'E[A_d]={EAdisc:.2f}={Adisc.mean():.2f}')
E[A_d]=3.39=3.39
The CDFs are plotted for the continuous interarrial time $A$ and the discretized interarrival time $A_d$.
t = np.linspace(0,4*EAdisc*deltaT,100)
plt.plot(t/deltaT, At(t), label=f'A ~ Exp($\lambda$) with $\lambda={lam:.2f}$')
plt.xlabel('interarrival time $A/\Delta t$')
plt.ylabel('CDF')
plt.grid()
plt.xlim([0,4*EAdisc])
Adisc.plotCDF(label='discretized $A_d$ ')
plt.legend()
plt.title(f'cont: $E[A]$={EA:.2f}, disc: $E[A_d]\cdot \Delta t$={EAdisc*deltaT:.2f}, $\Delta t={deltaT}$');
Note that the discreted exponential distribution yields a Geometric distribution with parameter $p=1-e^{-\lambda \Delta t}$.
$ \displaystyle A_d \sim \mathrm{GEOM}_1(1-e^{-\lambda \Delta t})$.
It is $ a_d(i+1) = A\big( (i+1) \Delta t\big) - A\big(i \Delta t\big) = e^{-\lambda i \Delta t} - e^{-\lambda (i+1) \Delta t} = e^{-\lambda i \Delta t}(1-e^{-\lambda \Delta t}) = (1-p)^i \cdot p $.
Remark: We could also arbitrarily define $a_d(i) = A\big( (i+1) \Delta t\big) - A\big(i \Delta t\big)$ yielding a GEOM$_0$ distribution with parameter $p$.
Adisc.plotPMF(label='$A_d$') # plot discretized interarrival time
D = GEOM(Adisc.mean(), m=1) # compare with GEOM1 distribution
D.plotPMF(label='GEOM1', linestyle='', marker='s')
# discretization of service time: same way
B = Bt(tdisc, EB=EB) # CDF of cont
B[-1] = 1
b = np.insert(np.diff(B), 0, 0)
Bdisc = DiscreteDistribution(np.arange(len(b)), b)
Bdisc.plotPMF(label='$B_d$') # plot discretized service time
plt.xlabel('time (s)')
plt.ylabel('probability')
plt.legend()
plt.xlim([-1, 10]);
plt.grid(which='major')
plt.title(title);
The discrete-time system function is the difference of the random variables: $ C = B-A $. The distribution is calculated by convolution.
$ \displaystyle c (k) = b (k) * a (-k) $
Cdisc = Bdisc - Adisc
Cdisc.plotPMF()
plt.xlabel('k')
plt.ylabel('c(k)=a(-k)*b(k)')
plt.xlim([-20,20])
plt.grid(which='major')
plt.title('system function c(k):'+title);
The waiting times of the (discretized) M/GI/1 are derived with the power method, until the difference of the waiting time distributions is below a threshold. For this purpose, the difference in the expected value of the waiting time in the iteration steps is computed and compared with a threshold value $ \epsilon $.
The expected value of the waiting time $ E [W_ {disc}] $ with the discrete-time analysis is compared with the exact value $ E [W] $ (for M / G / 1 waiting systems). In addition, Kingman's approximations for the expected value or an upper bound for the expected value are given.
Wn1 = DET(0) # empty system
Wn = DET(1) # just for initialization
# power method
while Wn != Wn1: # comparison based on means of the distributions
Wn = Wn1
Wn1 = max( Wn+Cdisc ,0)
Wdisc = Wn1
EW = EB * rho*(1+cB**2)/(2*(1-rho)) # exact formula for M/GI/1 delay system
#%% Kingman
Kingman = kingman(EA, cA, EB, cB) # Kingman approximation for mean waiting time
Upper = EA*(cA**2+rho**2*cB**2)/2/(1-rho) # upper bound of mean waiting time
print(f'Exact mean waiting time for M/GI/1 system: E[W]={EW:.4f}')
print(f'Kingman approximation: E[W_approx]={Kingman:.4f}')
print(f'Upper bound of the mean waiting time: {Upper:.4f}')
print(f'Discretized GI/GI/1 system: E[W_disc]*Delta t={Wdisc.mean()*deltaT:.4f}')
Exact mean waiting time for M/GI/1 system: E[W]=1.5556 Kingman approximation: E[W_approx]=1.5556 Upper bound of the mean waiting time: 2.7698 Discretized GI/GI/1 system: E[W_disc]*Delta t=1.4508
The cumulative distribution function of the waiting time (through the time-discrete analysis and discretization of the distributions) is now compared with the exact solution (M/GI/1 and Laplace transformation). For lower values of $ \Delta t $, the results become more precise at the expense of more computational effort.
Additionally, the exact CDF of the M/M/1 system is plotted. You may play with the parameter $k$; for $k=1$, the script computes M/M/1, for $k>1$, we have the M/E$_k$/1 queue.
t_lt = np.linspace(1e-9,10,100)
Wt = np.array([W(ti) for ti in t_lt])
plt.plot(t_lt/deltaT, Wt, label='M/GI/1 (numerical inverse Laplace)')
Wdisc.plotCDF(label='discretized M/GI/1')
# Exact solution for M/M/1: CDF
def Wmm1(t, lam=lam, EB=EB):
mu = 1/EB
rho = lam/mu
return 1-rho*np.exp(-mu*(1-rho)*t)
plt.plot(t_lt/deltaT, Wmm1(t_lt), label='M/M/1 (exact with CTMC)')
plt.xlabel('waiting time $t/\Delta t$')
plt.ylabel('CDF')
plt.grid(which='major')
plt.xlim([-1, 10])
plt.legend();
