I have created a simulation for an M/M/1 queuing system in Python using Simpy.
Simulations
The code of the simulation is just one class triggering two processes (in the Simpy sense of process):
- One for generating traffic.
- The other for describing the server.
import simpy
from random import expovariate
class QueueSys:
def __init__(self, env, lambd: float, mu: float):
self.env = env
self.__lambd = lambd # Traffic input rate (frequency)
self.__mu = mu # Server processing rate (frequency)
self.__counter = 0
self.__queue = []
self.__tracked_job_response_time = {} # Per each job, the times when they enter and leave the queue
self.__server = env.process(self.__serve())
self.__generator = env.process(self.__generate())
def get_stats(self):
result: list[float] = []
for job_id, job_times in self.__tracked_job_response_time.items():
t_start, t_end = job_times
if t_start > 0 and t_end > 0:
result.append(t_end - t_start)
return result
def __create_job(self):
job_name = f"J{self.__counter}"
self.__counter += 1
return job_name
def __enqueue(self, job: str):
self.__queue.append(job)
self.__track_job_start_time(job)
def __pop_job(self) -> str | None:
return self.__queue.pop(0) if len(self.__queue) > 0 else None
def __complete_job(self, job: str):
self.__track_job_end_time(job)
def __get_timeout_delay(self, rate: float) -> float:
return expovariate(rate)
def __track_job_start_time(self, job: str):
start_t, end_t = self.__tracked_job_response_time[job] if job in self.__tracked_job_response_time else (-1, -1)
if start_t > 0:
raise SystemError(f"Job '{job}' has already a tracked START time, unexpected")
start_t = self.env.now
self.__tracked_job_response_time[job] = (start_t, end_t)
def __track_job_end_time(self, job: str):
start_t, end_t = self.__tracked_job_response_time[job] if job in self.__tracked_job_response_time else (-1, -1)
if end_t > 0:
raise SystemError(f"Job '{job}' has already a tracked END time, unexpected")
end_t = self.env.now
self.__tracked_job_response_time[job] = (start_t, end_t)
def __serve(self):
while True:
processed_job = self.__pop_job()
yield self.env.timeout(self.__get_timeout_delay(self.__mu))
if processed_job is not None:
self.__complete_job(processed_job)
print("processed job %s at %d - qlen: %d" % (processed_job, self.env.now, len(self.__queue)))
def __generate(self):
while True:
new_job = self.__create_job()
self.__enqueue(new_job)
print("generated job '%s' at %d" % (new_job, self.env.now))
yield self.env.timeout(self.__get_timeout_delay(self.__lambd))
Running and configuring the simulation
I run the simulation:
env = simpy.Environment()
queue_mm1 = QueueSys(env, 0.5, 0.9)
env.run(until=5000)
results = queue_mm1.get_stats()
- Inter-arrivals are exponentially distributed with rate $\lambda$.
- Serving times are exponentially distributed with rate $\mu$.
I make sure that $\lambda < \mu$ so that the queue is stable.
The results
When I plot the results:
import matplotlib.pyplot as plt
plt.xlabel("Response time ranges")
plt.ylabel("Frequency")
plt.title("Response times histogram")
plt.hist(results, 300, density=True, facecolor='m')
plt.show()
I get this distribution:
The problem
For an M/M/1 queue, the response time of the system should be:
$$ E[T] = (\mu - \lambda)^{-1} $$
So, theoretically, the mean response time should be: $E[T] = \frac{1}{0.4} = 2.5$.
The issue here is that I am expecting response times to be exponentially distributed with mean $\mu - \lambda$. However I get what seems to be a Poisson distribution.
And the mean of response times does not match either:
mean(results)
And get approx: $3.56$.
What am I doing wrong in my simulation?
Additional notes
If I plot the number of items in the system (I do this by sampling the system size, queue + server, at a regular cadence), then I get the correct distribution: Geometrical.
Note: I did not update the code for this part.
The average number of items in the system is theoretically:
$$ E[N] = \frac{\rho}{1-\rho} $$
So, theoretically, I should have $1.25$ items on average in the system, and that is the case fr the simulation when I get the mean number of customers.
So, the simulations are getting right the number of items in the queue, but something is wrong with response times.
Resources
You can find a publicly available runnable Jupyter notebook with this simulation at this Google Colab page.
Indeed for an M/M/1 queue with $\lambda<\mu$, the steady state average delay is $\overline{W}=\frac{1}{\mu-\lambda}$ and the steady state delay distribution satisfies $P[W>x] = e^{-(\mu-\lambda)x}$ for all $x>0$.
If you have inter-arrival times $\{Y_k\}_{k=1}^{\infty}$ and service times $\{X_k\}_{k=1}^{\infty}$ then you can simulate a G/G/1 queue with delays $\{W_n\}_{n=1}^{\infty}$ given by $W_1=X_1$ and $$W_{n+1}= \max[W_n-Y_{n+1},0]+X_{n+1} \quad \forall n \in \{2, 3, 4, ...\}$$
For the special case of an M/M/1 queue with $\lambda=3, \mu=4$, it takes about $10^6$ arrivals to get a close match to theory, see Python simulation here: https://trinket.io/python3/9c74d2ebd6
Results for $\lambda=3,\mu=4$, $10^5$ iterations:
Empirical $E[W]$: 1.0293
Theory $E[W]$: 1.00
Empirical $P[W>1]$: 0.37761
Theory $P[W>1]$: 0.36787
Results for $\lambda=3, \mu=4, 10^6$ iterations:
Empirical $E[W]$: 1.003
Theory $E[W]$: 1.00
Empirical $P[W>1]$: 0.36788
Theory $P[W>1]$: 0.36787
Regarding what you wrote for Poisson distribution: If $Z$ is a random variable with a Poisson distribution, then $Z$ must take integer values with prob 1. So if your delays are always integers then you are not simulating correctly (since an M/M/1 queue has real-valued delays).