#### Previous topic

Priority Queues without Preemption

#### Next topic

Queueing Simulations by means of Messages

# Simulation of a Single Server Queueing Process¶

All code in one file

This code demonstrates how to build a simulator for the single server queue from first principles. Let $$a_i$$ be the time interval between customer $$i$$ and $$i+1$$, and $$b_i$$ be the amount of service of customer $$i$$. The numbers $$a_i, i = 1, 2, \ldots$$ are assumed to be independent identically distributed random variables with distribution function $$A$$, and the $$b_i, i = 1, 2, \ldots$$ are i.i.d. R.V. with distribution $$B$$. We assume that the first customer arrives at time $$t=0$$.

The second section shows how to implement the same simulation with knowledge propagation.

The last section considers a server subject to a D-policy

## The Single Server¶

The code below is not the fastest possible, but it is easy to spot the points where to speed it up.

import numpy as np
import scipy
import scipy.stats as stats
from pylab import plot, show

scipy.random.seed(3)
num = 35
labda = 1.
A = stats.expon(scale = 1./labda)
B = stats.uniform(scale = 1.5)


Let labda be the arrival rate of customers and num the number of customers that arrive during the simulation.

Now generate the interarrival times $$a_i$$ and service times $$b_i$$.

a = A.rvs(num-1)
b = B.rvs(num)


If $$A_i$$ is the arrival time of the $$i$$ th customer, then clearly, $$A_i = A_{i-1} + a_{i-1}$$. Using that the first customer arrives at $$t = 0$$:

def arrivalTimes(a):
At = np.zeros(len(a)+1)
At[1:] = np.cumsum(a)
return At


The waiting time for the $$i$$ th customer satisfies the following recursion:

$W_{i} = \max\{0, W_{i-1} - a_{i-1} + b_{i-1}, 0\}$
def waitingTimes(a,b, W0 = 0):
W = np.zeros_like(b)
W[0] = W0
for i in range(1,len(b)):
W[i] = max(W[i-1] - a[i-1] + b[i-1], 0)
return W


The sojourn time of the $$i$$ th customer, i.e., the total time this customer spends in the system, is his/her waiting time plus the service time, that is, $$S_i = W_i + b_i$$.

def sojournTimes(W, b):
S = W + b
return S


The departure time of a customer is easy: it is the arrival time plus the sojourn time. Therefore

def departureTimes(At, S):
return  At + S


Plotting the virtual waiting time requires some more attention. First, the virtual waiting time perceived by the $$i$$ th customer is the precisely the waiting time of the $$i$$ th customer. The virtual waiting time then immediately increases to $$W_i + b_i$$. Finally, observe that when $$W_i = 0$$, the virtual waiting time should be equal to zero at the departure moment of customer $$i$$.

def virtualWaitingTime(At, W, S, D):
V = []
for i in range(len(W)):
if i>0 and W[i] <= 0:
V.append([D[i-1],0])
V.append([At[i], W[i]])
V.append([At[i], S[i]])
V = np.array(V)
return V

At = arrivalTimes(a)
W = waitingTimes(a,b)
S = sojournTimes(W,b)
D = departureTimes(At,S)
V = virtualWaitingTime(At,W, S, D)
plot(V[:,0], V[:,1])
show()


## Using Trellis¶

Observe that in the code above all relations can be cast in the form of recursions. We exploit this idea here. To do so, we need the concept of a trellis, see <http://peak.telecommunity.com/DevCenter/Trellis#toc>. You may be interested in comparing this code to Priority queues with trellis.

from peak.events import trellis

class Job(trellis.Component):
trellis.attrs(
a = 0,  # interarrival time (to previous job)
b = 0,  # processing time
prev = None, # previous job in queue
)

trellis.maintain.attrs(
arTime = lambda self: self.prev.arTime + self.a if self.prev else 0,
waitTime = lambda self: max(self.prev.sojourn - self.a if self.prev else 0, 0),
sojourn = lambda self: self.waitTime + self.b,
depTime = lambda self: self.arTime +  self.sojourn,
)

def __repr__(self):
ret = str(self.a) + " "
ret += str(self.b) + " "
ret += str(self.arTime) + " "
ret += str(self.waitTime) + " "
ret += str(self.sojourn) + " "
ret += str(self.depTime) + " "
return ret

prev = Job(a=0,b=4)
sim = [prev]
for i in range(10):
nxt = Job(a=1, b=8, prev = prev)
sim.append(nxt)
prev = nxt

for j in sim:
print j

print "\n"

sim[3].b = 30
print sim[10]

0 4 0 0 4 4
1 8 1 3 11 12
1 8 2 10 18 20
1 8 3 17 25 28
1 8 4 24 32 36
1 8 5 31 39 44
1 8 6 38 46 52
1 8 7 45 53 60
1 8 8 52 60 68
1 8 9 59 67 76
1 8 10 66 74 84

1 8 10 88 96 106

Observe that now everything is implemented as a recursion. For instance, the job’s arrival time is sum of the arrival time of the its predecessor and the interarrival time. The waiting time is the sum of the sojourn time of its predecessor and its own process time, compare the relation

$W_{i} = \max\{0, W_{i-1} - a_{i-1} + b_{i-1}, 0\}$

Note that the sojourn time of the $$i$$ th job is equal to $$W_{i-1}+b_{i-1}$$. Etcetera. I actually find the code above much clearer.

## The D-Policy¶

We now simulate the waiting time process when the server is subjected to the, so called, D policy. This policy starts the server when the waiting time exceeds the threshold $$D$$, and stops the server when the waiting time hits zero again.

def waitingTimeDPolicy(a,b, D = 0, W0 = 0):
W = np.zeros_like(b)
W[0] = W0
On  = True if W0 > D else False
for i in range(1,len(b)):
if On:
W[i] = max(W[i-1] - a[i-1] + b[i-1], 0)
if W[i] <= 0:
On = False
else:
W[i] = W[i-1] + b[i-1]
if W[i] > D:
On = True
return W

At = arrivalTimes(a)
W = waitingTimeDPolicy(a,b, D = 5)
S = sojournTimes(W,b)
D = departureTimes(At,S)
V = virtualWaitingTime(At,W, S, D)
plot(V[:,0], V[:,1])
show()


In a simple cost model it costs $$c$$ per unit time per unit waiting time to keep customers waiting, and it costs $$K$$ to switch the server on. The following function computes the cost for a specific sample path.

c = 1.
K = 3
def computeCost(V):
S = 0. # total costs
T = 0. # total time
epoch, height = V[:,0], V[:,1]
for i in range(len(epoch)-1):
dx = epoch[i+1] - epoch[i]
dy = (height[i]+height[i+1])/2.
if dy == 0:
S += K
S += dx*dy
T += dx
return S, T

S, T = computeCost(V)
meanCost = S/T
print S, T, meanCost

133.675008897 31.7836442049 4.20577980407`