Simulating an M/M/1 queue

mm1-queue.py

def simulate_queue(arrival_rate, service_rate, simulation_time, time_step=0.1, seed=42):
    np.random.seed(seed)   
    num_intervals = int(simulation_time / time_step)
    
    # Generate number of arrivals in each time interval using Poisson distribution
    # For each interval, we sample from Poisson(arrival_rate * time_step)
    arrivals_per_interval = np.random.poisson(
        lam=arrival_rate * time_step, 
        size=num_intervals
    )
    
    # Convert to arrival times by distributing arrivals uniformly within each interval
    arrival_times = []
    for interval, num_arrivals in enumerate(arrivals_per_interval):
        if num_arrivals > 0:
            # For each arrival in this interval, place it uniformly within the interval
            interval_arrivals = interval * time_step + np.random.uniform(0, time_step, num_arrivals)
            arrival_times.extend(interval_arrivals)
    
    arrival_times = sorted(arrival_times)
    
    service_times = np.random.exponential(scale=1.0/service_rate, size=len(arrival_times))
    
    wait_times = []  # Wait time for each customer
    server_available_time = 0  # Time when server becomes available
    events = []  # List of (time, event_type) tuples: 1 for arrival, -1 for departure
    
    for i, arrival_time in enumerate(arrival_times):
        events.append((arrival_time, 1))
        
        # Calculate wait time for this customer
        if arrival_time >= server_available_time:
            wait_time = 0
            service_start_time = arrival_time
        else:
            wait_time = server_available_time - arrival_time
            service_start_time = server_available_time
            
        wait_times.append(wait_time)
        server_available_time = service_start_time + service_times[i]
        events.append((service_start_time + service_times[i], -1))
    
    events.sort()
    
    # Calculate queue length over time
    queue_length = 0
    queue_length_over_time = []
    time_points = []
    
    for time, event_type in events:
        if time <= simulation_time:
            queue_length += event_type  # +1 for arrival, -1 for departure
            queue_length_over_time.append(queue_length)
            time_points.append(time)
    
    return wait_times, queue_length_over_time, time_points