Warm-Up Period in SimPy: Avoiding Initialisation Bias

def run_simulation_with_warmup(warmup_time, collection_time):
    env = simpy.Environment()
    # ... setup

    # Phase 1: Warm-up (no statistics collection)
    env.run(until=warmup_time)

    # Reset statistics here
    stats.reset()

    # Phase 2: Data collection
    env.run(until=warmup_time + collection_time)

    return stats

def plot_warmup_detection(time_series_df, metric='queue_length'):
    """Visual inspection for warm-up period."""
    plt.figure(figsize=(12, 5))
    plt.plot(time_series_df['time'], time_series_df[metric])
    plt.xlabel('Time')
    plt.ylabel(metric)
    plt.title('Identify where the metric stabilises')
    plt.show()

def detect_warmup_moving_average(time_series, window_size=50, threshold=0.05):
    """Detect warm-up using moving average stability."""
    ma = pd.Series(time_series).rolling(window=window_size).mean()

    # Calculate relative change in moving average
    changes = ma.diff().abs() / ma

    # Find where changes fall below threshold
    stable_index = (changes < threshold).idxmax()

    return stable_index

def mser_warmup(data, max_truncation=0.5):
    """
    Find optimal truncation point using MSER-5.

    Truncates where mean squared error is minimised.
    """
    n = len(data)
    max_d = int(n * max_truncation)

    best_d = 0
    best_mser = float('inf')

    for d in range(0, max_d, 5):  # MSER-5 uses batches of 5
        remaining = data[d:]
        if len(remaining) < 10:
            break
        mser = np.var(remaining, ddof=1) / len(remaining)
        if mser < best_mser:
            best_mser = mser
            best_d = d

    return best_d

class Statistics:
    def __init__(self):
        self.wait_times = []
        self.collection_enabled = False

    def reset(self):
        self.wait_times = []
        self.collection_enabled = True

    def record_wait(self, wait):
        if self.collection_enabled:
            self.wait_times.append(wait)

# Usage
stats = Statistics()
env.run(until=warmup_time)
stats.reset()  # Start collecting
env.run(until=warmup_time + collection_time)

def customer(env, resource, stats, warmup_complete_event):
    arrival = env.now

    with resource.request() as req:
        yield req
        wait = env.now - arrival
        yield env.timeout(service_time())

        # Only record if warm-up is complete
        if warmup_complete_event.triggered:
            stats.record_wait(wait)

# Set flag after warm-up
def warmup_controller(env, warmup_time, warmup_event):
    yield env.timeout(warmup_time)
    warmup_event.succeed()
    print(f"Warm-up complete at {env.now}")

env = simpy.Environment()
warmup_event = env.event()
env.process(warmup_controller(env, 100, warmup_event))

def analyse_results(all_stats, warmup_time):
    """Filter out warm-up period after simulation."""
    return [
        record for record in all_stats
        if record['arrival_time'] >= warmup_time
    ]

def estimate_warmup_length(pilot_run_stats, target_stability=0.05):
    """
    Estimate required warm-up based on pilot run.

    target_stability: maximum coefficient of variation in rolling mean
    """
    means = pd.Series(pilot_run_stats).expanding().mean()
    cv = means.std() / means.mean()

    # Find where CV stabilises
    for i, val in enumerate(cv):
        if val < target_stability:
            return i

    return len(pilot_run_stats)  # Didn't stabilise

import simpy
import random
import numpy as np
import matplotlib.pyplot as plt

class WarmupAwareSimulation:
    def __init__(self, warmup_time, collection_time):
        self.warmup_time = warmup_time
        self.collection_time = collection_time
        self.wait_times = []
        self.all_wait_times = []  # Including warm-up for analysis
        self.warmup_complete = False

    def run(self, seed=42):
        random.seed(seed)
        env = simpy.Environment()
        server = simpy.Resource(env, capacity=2)

        env.process(self.arrivals(env, server))
        env.process(self.warmup_controller(env))

        env.run(until=self.warmup_time + self.collection_time)

        return self.wait_times

    def warmup_controller(self, env):
        yield env.timeout(self.warmup_time)
        self.warmup_complete = True
        print(f"Warm-up complete at time {env.now}")

    def arrivals(self, env, server):
        i = 0
        while True:
            yield env.timeout(random.expovariate(1/5))
            env.process(self.customer(env, f"C{i}", server))
            i += 1

    def customer(self, env, name, server):
        arrival = env.now
        with server.request() as req:
            yield req
            wait = env.now - arrival
            yield env.timeout(random.expovariate(1/4))

            self.all_wait_times.append({'time': arrival, 'wait': wait})
            if self.warmup_complete:
                self.wait_times.append(wait)

    def plot_warmup_analysis(self):
        times = [r['time'] for r in self.all_wait_times]
        waits = [r['wait'] for r in self.all_wait_times]

        plt.figure(figsize=(12, 5))
        plt.scatter(times, waits, alpha=0.3, s=5)
        plt.axvline(self.warmup_time, color='red', linestyle='--',
                    label='Warm-up end')
        plt.xlabel('Arrival Time')
        plt.ylabel('Wait Time')
        plt.title('Wait Times Over Simulation')
        plt.legend()
        plt.show()

# Run
sim = WarmupAwareSimulation(warmup_time=200, collection_time=1000)
results = sim.run()
sim.plot_warmup_analysis()

print(f"Customers in warm-up: {len(sim.all_wait_times) - len(sim.wait_times)}")
print(f"Customers analysed: {len(sim.wait_times)}")
print(f"Mean wait (steady-state): {np.mean(sim.wait_times):.2f}")