1. 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.
\inputminted[label=server.py]{python}{../code/event_stacks/server.py}
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.
\inputminted[label=servers.py]{python}{../code/event_stacks/servers.py}
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 \mintinline{python}{property}, as this allows us to access them the same way as, for example, the departure time.
The \mintinline{python}{__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.
\inputminted[label=job.py]{python}{../code/event_stacks/job.py}
The Queue class also subclassses list so that it can act as a heap and use Job.__lt__ to select the next job from the queue.
\inputminted[label=queues.py]{python}{../code/event_stacks/queues.py}
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.
\inputminted[label=event.py]{python}{../code/event_stacks/event.py}
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 arbitrary but still pop them sorted in time, you understand the core of discrete-event simulation.
\inputminted[label=events.py]{python}{../code/event_stacks/events.py}
The Statistics class is a subclass of list allowing us to append the jobs in the order in which they depart from the server.
\inputminted[label=stats.py]{python}{../code/event_stacks/stats.py}
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.
\inputminted[label=simulator.py]{python}{../code/event_stacks/simulator.py}
Testing the code comes next.
In \cref{sec:mm1}, 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 \cref{sec:simul-cont-time}. \inputminted[label=test-mmc.py]{python}{../code/event_stacks/test_mmc.py}
Here is the output for \(10^4\) jobs. The results are quite OK, but the differences in \(\E W\) and \(\E Q\) are still quite large.
exec(open("/home/nicky/vakken/qts/stochastic_or/code/event_stacks/test_mmc.py").read())
W: 0.0441, mm2.EW=0.0409
J: 0.2952, mm2.EJ=0.2909
Q: 0.1467, mm2.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 situation will we such numbers of jobs?
from override import run_with_overrides
script = (
"/home/nicky/vakken/qts/stochastic_or/code/event_stacks/test_mmc.py"
)
ns = {"num": 1000000}
print(run_with_overrides(script, ns))
W: 0.0410, mm2.EW=0.0409
J: 0.2915, mm2.EJ=0.2909
Q: 0.1225, mm2.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.
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 ChatGPT: `What is the advantage of a python deque over a heapq. Discuss the algorithmic efficiency, too.’
Consider a multiserver 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.
This work is licensed by the University of Groningen under a
Creative Commons Attribution-ShareAlike 4.0 International License.