3.5 Discrete-event Simulation
The queueing and inventory systems in the previous sections were based on simple recursions, and therefore straightforward to implement. When we need to analyze more difficult situations, we need discrete-event simulation. In this section we build a simulator for a \(G/G/c\) queue with different server speeds and test it. In the next section, we apply this to some non-trivial queueing situations. We will see the power of Python classes.
The first class is the Server which has an id and a service rate.
The method __lt__ enables us to express a preference for a server with a faster service rate.
from dataclasses import dataclass
@dataclass
class Server:
id: int
rate: float
def __lt__(self, other):
return self.rate > other.rate
The Servers class represents the pool of free servers.
We subclass Servers from list so that it can be used by the heapq algorithms to push and pop servers.
When a job leaves, Servers assigns the next job in queue to a server based on the priority rule stated in Server.__lt__.
A server that becomes idle after finishing a job is pushed to the server pool.
The convenience methods num_free and is_server_available improve readability compared to using len directly.
from heapq import heappop, heappush
class Servers(list):
def push(self, server):
heappush(self, server)
def pop(self):
return heappop(self)
def num_free(self):
return len(self)
def is_server_available(self):
return self.num_free() > 0
The Job class speaks mostly for itself.
The service time equals the job's load divided by the rate of the server that serves the job.1
Since servers may have different rates, we can only compute the service time after the job is assigned to a server.
The job keeps a reference to its assigned server.
The other attributes are used for the statistics computations at the end of the simulation.
We implement the sojourn time and waiting time as a property, as this allows us to access them the same way as, for example, the departure time.
The __lt__ method determines which job to select from the queue when a service can start.
As it is now, we select jobs in FIFO sequence.
The __hash__ method assigns each job a unique id, enabling its use in a set or as a key in a dict.
from dataclasses import dataclass
from server import Server
@dataclass
class Job:
id: int
arrival_time: float
load: float
service_time: float = 0
departure_time: float = 0
q_length_at_arrival: int = 0
q_length_at_departure: int = 0
free_servers: int = 0
server: Server = None
@property
def sojourn_time(self):
return self.departure_time - self.arrival_time
@property
def waiting_time(self):
return self.sojourn_time - self.service_time
def __str__(self):
return (
f"{self.id}, load={self.load:.2f}, service={self.service_time:.2f},"
f" arr={self.arrival_time:.2f}, dep={self.departure_time:.2f}"
)
def __hash__(self):
return self.id
def __lt__(self, other):
return self.arrival_time < other.arrival_time
The Queue class also subclasses list so that it can act as a heap and use Job.__lt__ to select the next job from the queue.
from heapq import heappop, heappush
class Queue(list):
def pop(self):
return heappop(self)
def push(self, job):
heappush(self, job)
def length(self):
return len(self)
def is_empty(self):
return self.length() == 0
An Event contains a time and a job; the __lt__ method specifies how to order events.2
Clearly, this must be by time.
The ArrivalEvent and DepartureEvent allow us to distinguish the relevant type of event.
from dataclasses import dataclass
from job import Job
@dataclass
class Event:
time: float
job: Job
def __lt__(self, other):
return self.time < other.time
class ArrivalEvent(Event):
pass
class DepartureEvent(Event):
pass
The Events class is a heap queue and keeps the events chronologically sorted.3
Once you understand that the Events class offers the functionality to remove and insert events arbitrarily but still pop them sorted in time, you understand the core of discrete-event simulation.
from heapq import heappop, heappush
class Events(list):
def push(self, event):
heappush(self, event)
def pop(self):
return heappop(self)
def is_empty(self):
return len(self) == 0
The Statistics class is a subclass of list allowing us to append the jobs in the order in which they depart from the server.
from job import Job
class Statistics(list[Job]):
def push(self, job: Job):
self.append(job)
def num_jobs(self):
return len(self)
def mean_waiting_time(self):
tot = sum(job.waiting_time for job in self)
return tot / self.num_jobs()
def mean_sojourn_time(self):
tot = sum(job.sojourn_time for job in self)
return tot / self.num_jobs()
def mean_queue_length(self):
tot = sum(job.q_length_at_arrival for job in self)
return tot / self.num_jobs()
def mean_servers_free(self):
tot = sum(job.free_servers for job in self)
return tot / self.num_jobs()
The Simulator is the final class 4
Study it carefully.
and it stores the events, queue, statistics, and servers.
The attribute now points to the current simulation time.
When starting the simulation, we provide a list of jobs to the simulator, which pushes an ArrivalEvent for each new job onto the event stack.
When a job service starts, the serve_job method asks for a free server from the server pool5
At a later stage, we might have different server pools. We therefore add this as an argument.
and assigns it to the job.
Then it computes the job service time and departure time, and pushes a DepartureEvent to the event stack to inform the simulator that a job will leave at that time.
The run method checks whether events remain; if so, it pops the next event from the event stack.
It then updates now to the event time and retrieves the associated job.
If the event is an ArrivalEvent, the job is served if a server is free; otherwise, it is queued.
In case of a DepartureEvent, the server assigned to the job becomes free and is returned to the server pool.
The departing job is pushed to the statistics tracker.
If there are still jobs in the queue, a new job is taken into service.
import events
import queues
import server
import servers
import stats
from event import ArrivalEvent, DepartureEvent
from job import Job
class Simulation:
def __init__(self):
self.events = events.Events()
self.queue = queues.Queue()
self.stats = stats.Statistics()
self.servers = servers.Servers()
self.now = 0
def initialize_jobs(self, jobs):
for job in jobs:
self.events.push(ArrivalEvent(job.arrival_time, job))
def initialize_servers(self, server_rates):
for i, rate in enumerate(server_rates):
self.servers.push(server.Server(id=i, rate=rate))
def serve_job(self, job: Job, server_pool=None):
"Serve job from a free server of the server pool."
if server_pool is None:
server_pool = self.servers
server = server_pool.pop()
job.free_servers = server_pool.num_free()
job.server = server
job.service_time = job.load / server.rate
job.departure_time = self.now + job.service_time
self.events.push(DepartureEvent(job.departure_time, job))
def handle_arrival(self, job):
job.q_length_at_arrival = self.queue.length()
if self.servers.is_server_available():
self.serve_job(job)
else:
self.queue.push(job)
def handle_departure(self, job):
self.servers.push(job.server)
self.stats.push(job)
if not self.queue.is_empty():
self.serve_job(self.queue.pop())
job.q_length_at_departure = self.queue.length()
def run(self):
while not self.events.is_empty():
event = self.events.pop()
self.now, job = event.time, event.job
if isinstance(event, ArrivalEvent):
self.handle_arrival(job)
elif isinstance(event, DepartureEvent):
self.handle_departure(job)
else:
raise ValueError("Unknown event")
Testing the code comes next.
In Section 5.1, we derive exact results for the expected waiting time of the \(M/M/c\) queue, so this model provides a good test.
The library mmc.py 6
See GitHub.
contains the code for the \(M/M/c\) queue.
For the generation of the data for the simulation, we follow the same procedure as in Section 3.3.
import numpy as np
import mmc
from job import Job
from simulator import Simulation
def run_mmc(num, labda=3, mu=4, server_rates=[1, 1]):
rng = np.random.default_rng(3)
X = rng.exponential(scale=1 / labda, size=num)
X[0] = 0
A = X.cumsum()
loads = rng.exponential(scale=1 / mu, size=num)
jobs = [Job(id=i, arrival_time=A[i], load=loads[i]) for i in range(1, num)]
sim = Simulation()
sim.initialize_jobs(jobs)
sim.initialize_servers(server_rates=server_rates)
sim.run()
c = len(server_rates)
mm = mmc.MMC(labda, mu, c=c)
print(f"W: {sim.stats.mean_waiting_time():0.4f}, {mm.EW=:.4f}")
print(f"J: {sim.stats.mean_sojourn_time():0.4f}, {mm.EJ=:.4f}")
print(f"Q: {sim.stats.mean_queue_length():0.4f}, {mm.EQ=:.4f}")
Here is the output for \(10^4\) jobs. The results are reasonable, but the differences in \(\E \W\) and \(\E \QQ\) are still quite large.
W: 0.0441, mm.EW=0.0409
J: 0.2952, mm.EJ=0.2909
Q: 0.1467, mm.EQ=0.1227
Let us repeat the experiment, but now with \(10^6\) jobs. The results show a remarkable improvement.7 This brings up a philosophical problem. If we need to simulate with one million jobs to see the time-average statistics appear, what is the significance of such time-averages? In how many practical situations will we see such numbers of jobs?
W: 0.0410, mm.EW=0.0409
J: 0.2915, mm.EJ=0.2909
Q: 0.1225, mm.EQ=0.1227
In the next section we demonstrate with a few examples how flexible this simulation environment actually is.
Claim: This part of the Simulator class handles arrivals correctly.
class Simulation:
def handle_arrival(self, job):
job.q_length_at_arrival = self.queue.length()
if self.servers.is_server_available():
self.queue.push(job)
else:
self.serve_job(job)
Solution
Solution, for real
False. Compare with the correct code.
Claim: This part of the simulator class handles the serving of jobs correctly.
class Simulation:
def serve_job(self, job):
server = self.servers.pop()
job.server = server
job.service_time = job.load / server.rate
job.departure_time = self.now + job.service_time
job.queue_length = self.queue.length()
job.free_servers = self.servers.num_free()
self.events.push(DepartureEvent(job.departure_time, job))
Solution
Solution, for real
True.
An interesting variation. Claim: the job records the number of free servers as seen just prior to arrival. This claim is false because the code first pops a server from the pool of free servers, and then sets the \pythoninline{job.free_servers} attribute.
Claim: This part of the Simulator class handles departures correctly.
class Simulation:
def handle_departure(self, job):
self.servers.push(job.server)
self.stats.push(job)
if not self.queue.is_empty():
self.serve_job(self.queue.pop())
job.q_length_at_departure = self.queue.length()
Solution
Solution, for real
True.
Claim: In discrete-event simulation, time advances by a fixed increment each step.
Solution
Solution, for real
False. In discrete-event simulation, time jumps from one event to the next. The time increments are variable and determined by the event schedule.
Claim: In a discrete-event simulation framework based on object-oriented programming, subclassing allows us to reuse the event loop and only override specific behavior such as arrival handling.
Solution
Solution, for real
True. The base class implements the general event loop and common logic, while subclasses override methods like handle_arrival to model specific queueing policies.
Consult the website of RealPython to gain a better understanding of classes, dataclasses, heapq, and deque.
Solution
Solution, for real
Read the pages on RealPython.
There are some advantages to subclassing Queue from deque.
What are these advantages?
Solution
Solution, for real
Give this string to an LLM: `What is the advantage of a python deque over a heapq. Discuss the algorithmic efficiency, too.'
Consider a multi-server queue with \(m\) servers. Suppose that at some time \(t\) it happens that \(\As(t) - D(t) < m\) but \(A(t) - D(t) > m\). How can this occur?
Solution
Solution, for real
In this case, there are servers idling while there are still customers in queue. If such events occur, the system is not work-conserving.