import heapq class Event: ''' Store the properties of one event in the Schedule class defined below. Each event has a time at which it needs to run, a function to call when running the event, along with the arguments and keyword arguments to pass to that function. ''' def __init__(self, timestamp, function, *args, **kwargs): self.timestamp = timestamp self.function = function self.args = args self.kwargs = kwargs def __lt__(self, other): ''' This overloads the less-than operator in Python. We need it so the priority queue knows how to compare two events. We want events with earlier (smaller) times to go first. ''' return self.timestamp < other.timestamp def run(self, schedule): ''' Run an event by calling the function with its arguments and keyword arguments. The first argument to any event function is always the schedule in which events are being tracked. The schedule object can be used to add new events to the priority queue. ''' self.function(schedule, *self.args, **self.kwargs) class Schedule: ''' Implement an event schedule using a priority queue. You can add events and run the next event. The `now` attribute contains the time at which the last event was run. ''' def __init__(self): self.now = 0 # Keep track of the current simulation time self.priority_queue = [] # The priority queue of events to run def add_event_at(self, timestamp, function, *args, **kwargs): # Add an event to the schedule at a particular point in time. heapq.heappush( self.priority_queue, Event(timestamp, function, *args, **kwargs)) def add_event_after(self, interval, function, *args, **kwargs): # Add an event to the schedule after a specified time interval. self.add_event_at(self.now + interval, function, *args, **kwargs) def next_event_time(self): return self.priority_queue[0].timestamp def run_next_event(self): # Get the next event from the priority queue and run it. event = heapq.heappop(self.priority_queue) self.now = event.timestamp event.run(self) def __repr__(self): return ( f'Schedule() at time {self.now} ' + f'with {len(self.priority_queue)} events in the queue') def print_events(self): print(repr(self)) for event in sorted(self.priority_queue): print(f' {event.timestamp}: {event.function.__name__}')

# M/D/1 queue import scipy.stats as sts class Queue: def __init__(self, service_rate): # Store the deterministic service time for an M/D/1 queue self.service_time = 1 / service_rate # We start with an empty queue and the server not busy self.people_in_queue = 0 self.people_being_served = 0 def add_customer(self, schedule): # Add the customer to the queue self.people_in_queue += 1 if self.people_being_served < 1: # This customer can be served immediately schedule.add_event_after(0, self.start_serving_customer) def start_serving_customer(self, schedule): # Move the customer from the queue to a server self.people_in_queue -= 1 self.people_being_served += 1 # Schedule when the server will be done with the customer schedule.add_event_after( self.service_time, self.finish_serving_customer) def finish_serving_customer(self, schedule): # Remove the customer from the server self.people_being_served -= 1 if self.people_in_queue > 0: # There are more people in the queue so serve the next customer schedule.add_event_after(0, self.start_serving_customer) class BusSystem: def __init__(self, arrival_rate, service_rate): self.queue = Queue(service_rate) self.arrival_distribution = sts.expon(scale=1/arrival_rate) def add_customer(self, schedule): # Add this customer to the queue self.queue.add_customer(schedule) # Schedule when to add another customer schedule.add_event_after( self.arrival_distribution.rvs(), self.add_customer) def run(self, schedule): # Schedule when the first customer arrives schedule.add_event_after( self.arrival_distribution.rvs(), self.add_customer) def run_simulation(arrival_rate, service_rate, run_until): schedule = Schedule() bus_system = BusSystem(arrival_rate, service_rate) bus_system.run(schedule) while schedule.next_event_time() < run_until: schedule.run_next_event() return bus_system

# Run a short test bus_system = run_simulation(arrival_rate=0.8, service_rate=1, run_until=100) print(f'There are {bus_system.queue.people_in_queue} people in the queue')

# Iterating the simulation (100 trials) import matplotlib.pyplot as plt storage_100 = [] iteration = 100 for i in range(iteration): temp = run_simulation(arrival_rate=0.8, service_rate=1, run_until=100) storage_100.append(temp.queue.people_in_queue)

# Plotting the historgram plt.hist(storage_100, bins = 10) plt.xlabel("Queue length") plt.ylabel("Count")

# Confidence interval calculation import numpy as np data = storage_100 m = np.mean(data) print('Sample mean:', m) t = sts.sem(data) print('Standard error of the mean:', t) print('95% confidence interval of population mean:', [m + 1.96*t, m - 1.96*t])

# Iterating the simulation (1000 trials) import matplotlib.pyplot as plt storage_1000 = [] iteration = 1000 for i in range(iteration): temp = run_simulation(arrival_rate=0.8, service_rate=1, run_until=100) storage_1000.append(temp.queue.people_in_queue)

# Plotting the historgram plt.hist(storage_1000, bins = 10) plt.xlabel("Queue length") plt.ylabel("Count")

# Confidence interval calculation import numpy as np data = storage_1000 m = np.mean(data) print('Sample mean:', m) t = sts.sem(data) print('Standard error of the mean:', t) print('95% confidence interval of population mean:', [m + 1.96*t, m - 1.96*t])

# Iterating the simulation (1000 trials) import matplotlib.pyplot as plt storage_1 = [] storage_2 = [] storage_3 = [] storage_4 = [] storage_5 = [] storage_6 = [] storage_7 = [] iteration = 1000 # Run the simulations for i in range(iteration): temp_1 = run_simulation(arrival_rate=0.2, service_rate=1, run_until=100) storage_1.append(temp_1.queue.people_in_queue) temp_2 = run_simulation(arrival_rate=0.4, service_rate=1, run_until=100) storage_2.append(temp_2.queue.people_in_queue) temp_3 = run_simulation(arrival_rate=0.6, service_rate=1, run_until=100) storage_3.append(temp_3.queue.people_in_queue) temp_4 = run_simulation(arrival_rate=0.8, service_rate=1, run_until=100) storage_4.append(temp_4.queue.people_in_queue) temp_5 = run_simulation(arrival_rate=1, service_rate=1, run_until=100) storage_5.append(temp_5.queue.people_in_queue) temp_6 = run_simulation(arrival_rate=1.2, service_rate=1, run_until=100) storage_6.append(temp_6.queue.people_in_queue) temp_7 = run_simulation(arrival_rate=1.4, service_rate=1, run_until=100) storage_7.append(temp_7.queue.people_in_queue)

# Variables to store error plot x = [i/10 for i in range(2, 16, 2)] y = [] y_error = [] # Find the average queue length and error for each model m_1 = np.mean(storage_1) t_1 = sts.sem(storage_1) y.append(m_1) y_error.append(m_1 + 1.96*t_1) m_2 = np.mean(storage_2) t_2 = sts.sem(storage_2) y.append(m_2) y_error.append(m_2 + 1.96*t_2) m_3 = np.mean(storage_3) t_3 = sts.sem(storage_3) y.append(m_3) y_error.append(m_3 + 1.96*t_3) m_4 = np.mean(storage_4) t_4 = sts.sem(storage_4) y.append(m_4) y_error.append(m_4 + 1.96*t_4) m_5 = np.mean(storage_5) t_5 = sts.sem(storage_5) y.append(m_5) y_error.append(m_5 + 1.96*t_5) m_6 = np.mean(storage_6) t_6 = sts.sem(storage_6) y.append(m_6) y_error.append(m_6 + 1.96*t_6) m_7 = np.mean(storage_7) t_7 = sts.sem(storage_7) y.append(m_7) y_error.append(m_7 + 1.96*t_7)

# Error plot of average queue length and its 95% confidence interval plt.errorbar(x, y, y_error, color='black', marker='o', capsize=5, linestyle='--', linewidth=1, label="simulations") plt.xlabel("Arrival rate") plt.ylabel("Average queue length") plt.plot([0, 1.4], [5, 5], color='red', linestyle='-', linewidth=1, label = "queue length = 5") plt.legend() plt.show()