Industry Example SimPy Simulations - School of Simulation

Steel Production Plant Simulation

This simulation models a steel production facility with multiple furnaces processing steel batches. The simulation demonstrates:

Resource allocation and contention (furnaces)
Process times with variability
Queue management and waiting times
Production throughput analysis

Adjust the parameters below to see how different configurations affect production efficiency and resource utilization.

Simulation Parameters

Number of Furnaces

Number of Batches

Mean Process Time (hours)

Mean Arrival Interval (hours)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_FURNACES = {NUM_FURNACES}
NUM_BATCHES = {NUM_BATCHES}
PROCESS_TIME_MEAN = {PROCESS_TIME_MEAN}
ARRIVAL_INTERVAL_MEAN = {ARRIVAL_INTERVAL_MEAN}
SIM_TIME = 100

# Data collection
results = []

def steel_batch(env, name, furnaces):
    """Process a steel batch through the production line."""
    arrival_time = env.now
    
    # Request a furnace
    with furnaces.request() as request:
        yield request
        wait_time = env.now - arrival_time
        
        # Process the steel
        process_time = random.expovariate(1.0 / PROCESS_TIME_MEAN)
        yield env.timeout(process_time)
        
        # Record results
        results.append({
            'batch': name,
            'arrival_time': arrival_time,
            'wait_time': wait_time,
            'process_time': process_time,
            'completion_time': env.now
        })
        
        print(f"Batch {name}: Arrived at {arrival_time:.2f}, "
              f"waited {wait_time:.2f}, processed in {process_time:.2f}")

def batch_generator(env, furnaces):
    """Generate steel batches at random intervals."""
    for i in range(NUM_BATCHES):
        yield env.timeout(random.expovariate(1.0 / ARRIVAL_INTERVAL_MEAN))
        env.process(steel_batch(env, i, furnaces))

# Run simulation
env = simpy.Environment()
furnaces = simpy.Resource(env, capacity=NUM_FURNACES)
env.process(batch_generator(env, furnaces))
env.run(until=SIM_TIME)

# Analysis
df = pd.DataFrame(results)
print(f"\n=== Simulation Results ===")
print(f"Total batches processed: {len(df)}")
print(f"Average wait time: {df['wait_time'].mean():.2f} hours")
print(f"Average process time: {df['process_time'].mean():.2f} hours")
print(f"Average total time: {(df['wait_time'] + df['process_time']).mean():.2f} hours")
print(f"Furnace utilization: {(df['process_time'].sum() / (NUM_FURNACES * SIM_TIME) * 100):.1f}%")

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Wait times over time
axes[0, 0].plot(df['arrival_time'], df['wait_time'], 'o', alpha=0.7)
axes[0, 0].set_xlabel('Arrival Time (hours)')
axes[0, 0].set_ylabel('Wait Time (hours)')
axes[0, 0].set_title('Wait Times Throughout Simulation')
axes[0, 0].grid(True, alpha=0.3)

# 2. Process time distribution
axes[0, 1].hist(df['process_time'], bins=20, alpha=0.7, color='green', edgecolor='black')
axes[0, 1].set_xlabel('Process Time (hours)')
axes[0, 1].set_ylabel('Frequency')
axes[0, 1].set_title('Process Time Distribution')
axes[0, 1].axvline(df['process_time'].mean(), color='red', linestyle='--', label=f'Mean: {df["process_time"].mean():.2f}')
axes[0, 1].legend()

# 3. Queue length over time (approximation)
queue_events = []
for _, row in df.iterrows():
    queue_events.append((row['arrival_time'], 1))  # Arrival
    queue_events.append((row['arrival_time'] + row['wait_time'], -1))  # Start processing
queue_events.sort()

times, queue_length = [], []
current_queue = 0
for time, change in queue_events:
    times.append(time)
    current_queue += change
    queue_length.append(current_queue)

axes[1, 0].step(times, queue_length, where='post', alpha=0.7)
axes[1, 0].set_xlabel('Time (hours)')
axes[1, 0].set_ylabel('Queue Length')
axes[1, 0].set_title('Queue Length Over Time')
axes[1, 0].grid(True, alpha=0.3)

# 4. Throughput analysis
hourly_completions = df.groupby(df['completion_time'] // 1).size()
axes[1, 1].bar(hourly_completions.index, hourly_completions.values, alpha=0.7)
axes[1, 1].set_xlabel('Hour')
axes[1, 1].set_ylabel('Batches Completed')
axes[1, 1].set_title('Hourly Production Throughput')
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")  # Special marker for plot data

Output

Click "Run Simulation" to see results...

Shipping Logistics Simulation

This simulation models a shipping port with multiple berths handling container ships. The simulation demonstrates:

Port berth utilization
Ship arrival patterns and queuing
Loading/unloading time variability
Port throughput optimization

Experiment with different numbers of berths and ship arrival rates to optimize port operations.

Simulation Parameters

Number of Berths

Number of Ships

Mean Unload Time (hours)

Mean Arrival Interval (hours)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_BERTHS = {NUM_BERTHS}
NUM_SHIPS = {NUM_SHIPS}
UNLOAD_TIME_MEAN = {UNLOAD_TIME_MEAN}
ARRIVAL_INTERVAL_MEAN = {ARRIVAL_INTERVAL_MEAN}
SIM_TIME = 200

# Data collection
ship_data = []
berth_utilization = {i: [] for i in range(NUM_BERTHS)}

def ship_process(env, ship_id, berths):
    """Process a ship through the port."""
    arrival_time = env.now
    
    # Request a berth
    with berths.request() as request:
        yield request
        wait_time = env.now - arrival_time
        berth_id = berths.users.index(request)
        
        # Unload the ship
        unload_time = random.expovariate(1.0 / UNLOAD_TIME_MEAN)
        start_unload = env.now
        yield env.timeout(unload_time)
        
        # Record berth utilization
        berth_utilization[berth_id].append({
            'start': start_unload,
            'end': env.now,
            'ship': ship_id
        })
        
        # Record ship data
        ship_data.append({
            'ship_id': ship_id,
            'arrival_time': arrival_time,
            'wait_time': wait_time,
            'unload_time': unload_time,
            'departure_time': env.now,
            'berth_used': berth_id
        })
        
        print(f"Ship {ship_id}: Arrived at {arrival_time:.2f}, "
              f"waited {wait_time:.2f}, unloaded in {unload_time:.2f} at berth {berth_id}")

def ship_generator(env, berths):
    """Generate ships arriving at the port."""
    for i in range(NUM_SHIPS):
        yield env.timeout(random.expovariate(1.0 / ARRIVAL_INTERVAL_MEAN))
        env.process(ship_process(env, i, berths))

# Run simulation
env = simpy.Environment()
berths = simpy.Resource(env, capacity=NUM_BERTHS)
env.process(ship_generator(env, berths))
env.run(until=SIM_TIME)

# Analysis
df = pd.DataFrame(ship_data)
print(f"\n=== Port Performance Metrics ===")
print(f"Total ships processed: {len(df)}")
print(f"Average wait time: {df['wait_time'].mean():.2f} hours")
print(f"Average unload time: {df['unload_time'].mean():.2f} hours")
print(f"Average turnaround time: {(df['wait_time'] + df['unload_time']).mean():.2f} hours")

# Calculate berth utilization
total_berth_busy_time = sum(df['unload_time'])
berth_utilization_pct = (total_berth_busy_time / (NUM_BERTHS * SIM_TIME) * 100)
print(f"Overall berth utilization: {berth_utilization_pct:.1f}%")

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Ship wait times over time
axes[0, 0].scatter(df['arrival_time'], df['wait_time'], alpha=0.6, s=50)
axes[0, 0].set_xlabel('Arrival Time (hours)')
axes[0, 0].set_ylabel('Wait Time (hours)')
axes[0, 0].set_title('Ship Wait Times Throughout Simulation')
axes[0, 0].grid(True, alpha=0.3)

# 2. Berth utilization timeline
colors = plt.cm.tab10(np.linspace(0, 1, NUM_BERTHS))
for berth_id, periods in berth_utilization.items():
    for period in periods:
        axes[0, 1].barh(berth_id, period['end'] - period['start'], 
                       left=period['start'], height=0.8, 
                       color=colors[berth_id], alpha=0.7)
axes[0, 1].set_xlabel('Time (hours)')
axes[0, 1].set_ylabel('Berth ID')
axes[0, 1].set_title('Berth Utilization Timeline')
axes[0, 1].set_ylim(-0.5, NUM_BERTHS - 0.5)
axes[0, 1].grid(True, alpha=0.3, axis='x')

# 3. Turnaround time distribution
turnaround_times = df['wait_time'] + df['unload_time']
axes[1, 0].hist(turnaround_times, bins=20, alpha=0.7, color='coral', edgecolor='black')
axes[1, 0].set_xlabel('Turnaround Time (hours)')
axes[1, 0].set_ylabel('Frequency')
axes[1, 0].set_title('Ship Turnaround Time Distribution')
axes[1, 0].axvline(turnaround_times.mean(), color='red', linestyle='--', 
                   label=f'Mean: {turnaround_times.mean():.2f}')
axes[1, 0].legend()

# 4. Queue length analysis
queue_events = []
for _, row in df.iterrows():
    queue_events.append((row['arrival_time'], 1))
    queue_events.append((row['arrival_time'] + row['wait_time'], -1))
queue_events.sort()

times, queue_length = [0], [0]
current_queue = 0
for time, change in queue_events:
    if time <= SIM_TIME:
        times.append(time)
        current_queue += change
        queue_length.append(current_queue)

axes[1, 1].step(times, queue_length, where='post', alpha=0.7, linewidth=2)
axes[1, 1].fill_between(times, queue_length, step='post', alpha=0.3)
axes[1, 1].set_xlabel('Time (hours)')
axes[1, 1].set_ylabel('Ships in Queue')
axes[1, 1].set_title('Port Queue Length Over Time')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")  # Special marker for plot data

Output

Click "Run Simulation" to see results...

Satellite Communication Network

This simulation models a satellite communication network managing data transmission requests. The simulation demonstrates:

Satellite channel allocation and bandwidth management
Data packet transmission with varying priorities
Network congestion and latency analysis
Ground station coordination

Experiment with different numbers of satellites and data loads to optimize network performance.

Simulation Parameters

Number of Satellites

Number of Transmission Requests

Mean Transmission Time (minutes)

Mean Request Interval (minutes)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_SATELLITES = {NUM_SATELLITES}
NUM_REQUESTS = {NUM_REQUESTS}
TRANSMISSION_TIME_MEAN = {TRANSMISSION_TIME_MEAN}
REQUEST_INTERVAL_MEAN = {REQUEST_INTERVAL_MEAN}
SIM_TIME = 300  # 5 hours

# Data collection
transmission_data = []
satellite_utilization = {i: [] for i in range(NUM_SATELLITES)}

def data_transmission(env, request_id, satellites, priority):
    """Process a data transmission request through the satellite network."""
    arrival_time = env.now
    
    # Request a satellite channel
    with satellites.request(priority=priority) as request:
        yield request
        wait_time = env.now - arrival_time
        satellite_id = satellites.users.index(request) % NUM_SATELLITES
        
        # Transmit data
        transmission_time = random.expovariate(1.0 / TRANSMISSION_TIME_MEAN)
        start_transmission = env.now
        yield env.timeout(transmission_time)
        
        # Record satellite utilization
        satellite_utilization[satellite_id].append({
            'start': start_transmission,
            'end': env.now,
            'request': request_id
        })
        
        # Record transmission data
        transmission_data.append({
            'request_id': request_id,
            'priority': priority,
            'arrival_time': arrival_time,
            'wait_time': wait_time,
            'transmission_time': transmission_time,
            'completion_time': env.now,
            'satellite_used': satellite_id,
            'data_volume': random.uniform(10, 100)  # MB
        })
        
        priority_name = 'High' if priority == 1 else 'Medium' if priority == 2 else 'Low'
        print(f"Request {request_id} ({priority_name} priority): Arrived at {arrival_time:.2f}, "
              f"waited {wait_time:.2f}, transmitted in {transmission_time:.2f} via Satellite {satellite_id}")

def request_generator(env, satellites):
    """Generate data transmission requests with varying priorities."""
    for i in range(NUM_REQUESTS):
        yield env.timeout(random.expovariate(1.0 / REQUEST_INTERVAL_MEAN))
        # Assign priority: 1 (high), 2 (medium), 3 (low)
        priority = random.choices([1, 2, 3], weights=[0.2, 0.5, 0.3])[0]
        env.process(data_transmission(env, i, satellites, priority))

# Run simulation
env = simpy.Environment()
satellites = simpy.PriorityResource(env, capacity=NUM_SATELLITES)
env.process(request_generator(env, satellites))
env.run(until=SIM_TIME)

# Analysis
df = pd.DataFrame(transmission_data)
print(f"\n=== Satellite Network Performance ===")
print(f"Total requests processed: {len(df)}")
print(f"Average wait time: {df['wait_time'].mean():.2f} minutes")
print(f"Average transmission time: {df['transmission_time'].mean():.2f} minutes")
print(f"Average turnaround time: {(df['wait_time'] + df['transmission_time']).mean():.2f} minutes")

# Priority analysis
for priority in [1, 2, 3]:
    priority_df = df[df['priority'] == priority]
    if len(priority_df) > 0:
        priority_name = 'High' if priority == 1 else 'Medium' if priority == 2 else 'Low'
        print(f"\n{priority_name} Priority - Avg wait: {priority_df['wait_time'].mean():.2f} min")

# Calculate satellite utilization
total_busy_time = sum(df['transmission_time'])
utilization_pct = (total_busy_time / (NUM_SATELLITES * SIM_TIME) * 100)
print(f"\nOverall satellite utilization: {utilization_pct:.1f}%")

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Wait times by priority
priority_names = {1: 'High', 2: 'Medium', 3: 'Low'}
df['priority_name'] = df['priority'].map(priority_names)
sns.boxplot(data=df, x='priority_name', y='wait_time', ax=axes[0, 0])
axes[0, 0].set_xlabel('Priority Level')
axes[0, 0].set_ylabel('Wait Time (minutes)')
axes[0, 0].set_title('Wait Times by Priority Level')

# 2. Satellite utilization timeline
colors = plt.cm.viridis(np.linspace(0, 1, NUM_SATELLITES))
for sat_id, periods in satellite_utilization.items():
    for period in periods:
        axes[0, 1].barh(sat_id, period['end'] - period['start'], 
                       left=period['start'], height=0.8, 
                       color=colors[sat_id], alpha=0.7)
axes[0, 1].set_xlabel('Time (minutes)')
axes[0, 1].set_ylabel('Satellite ID')
axes[0, 1].set_title('Satellite Channel Utilization Timeline')
axes[0, 1].set_ylim(-0.5, NUM_SATELLITES - 0.5)
axes[0, 1].grid(True, alpha=0.3, axis='x')

# 3. Data throughput over time
df['hour'] = df['completion_time'] // 60
hourly_data = df.groupby('hour')['data_volume'].sum()
axes[1, 0].bar(hourly_data.index, hourly_data.values, alpha=0.7, color='teal')
axes[1, 0].set_xlabel('Hour')
axes[1, 0].set_ylabel('Data Volume (MB)')
axes[1, 0].set_title('Hourly Data Throughput')
axes[1, 0].grid(True, alpha=0.3, axis='y')

# 4. Network load over time
load_events = []
for _, row in df.iterrows():
    load_events.append((row['arrival_time'], 1))
    load_events.append((row['completion_time'], -1))
load_events.sort()

times, active_transmissions = [0], [0]
current_load = 0
for time, change in load_events:
    if time <= SIM_TIME:
        times.append(time)
        current_load += change
        active_transmissions.append(current_load)

axes[1, 1].plot(times, active_transmissions, alpha=0.7, linewidth=2)
axes[1, 1].fill_between(times, active_transmissions, alpha=0.3)
axes[1, 1].axhline(y=NUM_SATELLITES, color='red', linestyle='--', label=f'Capacity: {NUM_SATELLITES}')
axes[1, 1].set_xlabel('Time (minutes)')
axes[1, 1].set_ylabel('Active Transmissions')
axes[1, 1].set_title('Network Load Over Time')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")  # Special marker for plot data

Output

Click "Run Simulation" to see results...

Aircraft Turnaround Operations

This simulation models airport ground operations for aircraft turnaround between flights. The simulation demonstrates:

Gate allocation and aircraft parking
Parallel ground services (fueling, catering, cleaning, baggage)
Turnaround time optimization
Resource scheduling and coordination

Adjust parameters to minimize turnaround times while managing limited ground resources.

Simulation Parameters

Number of Gates

Number of Flights

Target Turnaround Time (minutes)

Mean Arrival Interval (minutes)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_GATES = {NUM_GATES}
NUM_FLIGHTS = {NUM_FLIGHTS}
TARGET_TURNAROUND = {TARGET_TURNAROUND}
ARRIVAL_INTERVAL_MEAN = {ARRIVAL_INTERVAL_MEAN}
SIM_TIME = 1440  # 24 hours

# Ground service resources (shared across all gates)
NUM_FUEL_TRUCKS = 4
NUM_CATERING_TRUCKS = 3
NUM_CLEANING_CREWS = 5
NUM_BAGGAGE_CREWS = 4

# Service time parameters (in minutes)
FUELING_TIME = 20
CATERING_TIME = 15
CLEANING_TIME = 25
BAGGAGE_TIME = 30

# Data collection
flight_data = []
gate_utilization = {i: [] for i in range(NUM_GATES)}

def ground_service(env, name, duration, resource, flight_id):
    """Perform a ground service task."""
    with resource.request() as request:
        yield request
        start_time = env.now
        yield env.timeout(duration + random.uniform(-5, 5))  # Add variability
        return env.now - start_time

def aircraft_turnaround(env, flight_id, gates, fuel_trucks, catering_trucks, 
                       cleaning_crews, baggage_crews):
    """Process aircraft turnaround operations."""
    arrival_time = env.now
    
    # Request a gate
    with gates.request() as gate_request:
        yield gate_request
        gate_wait = env.now - arrival_time
        gate_id = gates.users.index(gate_request) % NUM_GATES
        turnaround_start = env.now
        
        # Perform ground services in parallel
        services = [
            env.process(ground_service(env, 'Fueling', FUELING_TIME, fuel_trucks, flight_id)),
            env.process(ground_service(env, 'Catering', CATERING_TIME, catering_trucks, flight_id)),
            env.process(ground_service(env, 'Cleaning', CLEANING_TIME, cleaning_crews, flight_id)),
            env.process(ground_service(env, 'Baggage', BAGGAGE_TIME, baggage_crews, flight_id))
        ]
        
        # Wait for all services to complete
        service_times = yield env.all_of(services)
        
        turnaround_time = env.now - turnaround_start
        total_time = env.now - arrival_time
        
        # Record gate utilization
        gate_utilization[gate_id].append({
            'start': turnaround_start,
            'end': env.now,
            'flight': flight_id
        })
        
        # Record flight data
        flight_data.append({
            'flight_id': flight_id,
            'arrival_time': arrival_time,
            'gate_wait': gate_wait,
            'turnaround_time': turnaround_time,
            'total_time': total_time,
            'departure_time': env.now,
            'gate_used': gate_id,
            'on_time': turnaround_time <= TARGET_TURNAROUND
        })
        
        status = '✓ On-time' if turnaround_time <= TARGET_TURNAROUND else '✗ Delayed'
        print(f"Flight {flight_id:03d}: Gate {gate_id}, Turnaround {turnaround_time:.1f} min {status}")

def flight_generator(env, gates, fuel_trucks, catering_trucks, cleaning_crews, baggage_crews):
    """Generate arriving flights."""
    for i in range(NUM_FLIGHTS):
        yield env.timeout(random.expovariate(1.0 / ARRIVAL_INTERVAL_MEAN))
        env.process(aircraft_turnaround(env, i, gates, fuel_trucks, catering_trucks,
                                       cleaning_crews, baggage_crews))

# Run simulation
env = simpy.Environment()
gates = simpy.Resource(env, capacity=NUM_GATES)
fuel_trucks = simpy.Resource(env, capacity=NUM_FUEL_TRUCKS)
catering_trucks = simpy.Resource(env, capacity=NUM_CATERING_TRUCKS)
cleaning_crews = simpy.Resource(env, capacity=NUM_CLEANING_CREWS)
baggage_crews = simpy.Resource(env, capacity=NUM_BAGGAGE_CREWS)

env.process(flight_generator(env, gates, fuel_trucks, catering_trucks, 
                           cleaning_crews, baggage_crews))
env.run(until=SIM_TIME)

# Analysis
df = pd.DataFrame(flight_data)
print(f"\n=== Airport Operations Performance ===")
print(f"Total flights processed: {len(df)}")
print(f"Average gate wait: {df['gate_wait'].mean():.2f} minutes")
print(f"Average turnaround time: {df['turnaround_time'].mean():.2f} minutes")
print(f"On-time performance: {(df['on_time'].sum() / len(df) * 100):.1f}%")
print(f"Flights meeting target: {df['on_time'].sum()} / {len(df)}")

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Turnaround time distribution
axes[0, 0].hist(df['turnaround_time'], bins=20, alpha=0.7, color='skyblue', edgecolor='black')
axes[0, 0].axvline(TARGET_TURNAROUND, color='red', linestyle='--', 
                   label=f'Target: {TARGET_TURNAROUND} min')
axes[0, 0].axvline(df['turnaround_time'].mean(), color='green', linestyle='--', 
                   label=f'Mean: {df["turnaround_time"].mean():.1f} min')
axes[0, 0].set_xlabel('Turnaround Time (minutes)')
axes[0, 0].set_ylabel('Number of Flights')
axes[0, 0].set_title('Aircraft Turnaround Time Distribution')
axes[0, 0].legend()

# 2. Gate utilization Gantt chart
colors = plt.cm.Set3(np.linspace(0, 1, NUM_GATES))
for gate_id, periods in gate_utilization.items():
    for period in periods:
        axes[0, 1].barh(gate_id, period['end'] - period['start'], 
                       left=period['start'], height=0.8, 
                       color=colors[gate_id], alpha=0.7)
axes[0, 1].set_xlabel('Time (minutes)')
axes[0, 1].set_ylabel('Gate Number')
axes[0, 1].set_title('Gate Utilization Schedule')
axes[0, 1].set_ylim(-0.5, NUM_GATES - 0.5)
axes[0, 1].grid(True, alpha=0.3, axis='x')

# 3. On-time performance over time
df['hour'] = df['arrival_time'] // 60
hourly_performance = df.groupby('hour')['on_time'].mean() * 100
axes[1, 0].plot(hourly_performance.index, hourly_performance.values, 
                marker='o', linewidth=2, markersize=8)
axes[1, 0].axhline(y=80, color='green', linestyle='--', alpha=0.5, label='80% target')
axes[1, 0].set_xlabel('Hour of Day')
axes[1, 0].set_ylabel('On-time Performance (%)')
axes[1, 0].set_title('On-time Performance by Hour')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].legend()
axes[1, 0].set_ylim(0, 105)

# 4. Resource utilization overview
resources = ['Gates', 'Fuel Trucks', 'Catering', 'Cleaning', 'Baggage']
utilizations = [
    (df['turnaround_time'].sum() / (NUM_GATES * SIM_TIME) * 100),
    60,  # Estimated
    50,  # Estimated
    70,  # Estimated
    65   # Estimated
]
bars = axes[1, 1].bar(resources, utilizations, alpha=0.7, 
                      color=['blue', 'orange', 'green', 'red', 'purple'])
axes[1, 1].axhline(y=80, color='red', linestyle='--', alpha=0.5)
axes[1, 1].set_ylabel('Utilization (%)')
axes[1, 1].set_title('Ground Resource Utilization')
axes[1, 1].set_ylim(0, 100)
axes[1, 1].grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar in bars:
    height = bar.get_height()
    axes[1, 1].text(bar.get_x() + bar.get_width()/2., height,
                    f'{height:.0f}%', ha='center', va='bottom')

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")  # Special marker for plot data

Output

Click "Run Simulation" to see results...

Military Supply Chain

This simulation models a military supply chain network with multiple forward operating bases (FOBs). The simulation demonstrates:

Multi-echelon inventory management
Supply convoy scheduling and routing
Emergency resupply prioritization
Resource allocation under constraints

Experiment with different supply strategies to maintain operational readiness across all bases.

Simulation Parameters

Number of FOBs

Number of Supply Convoys

Daily Consumption Rate (%)

Mean Delivery Time (days)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_FOBS = {NUM_FOBS}
NUM_CONVOYS = {NUM_CONVOYS}
DAILY_CONSUMPTION_RATE = {DAILY_CONSUMPTION_RATE} / 100
DELIVERY_TIME_MEAN = {DELIVERY_TIME_MEAN} * 24  # Convert to hours
SIM_TIME = 720  # 30 days in hours

# Supply categories
SUPPLY_TYPES = ['Ammunition', 'Fuel', 'Food', 'Medical']

# Data collection
supply_data = []
convoy_missions = []
base_status = {i: {'inventory': {s: 100 for s in SUPPLY_TYPES}} for i in range(NUM_FOBS)}

class FOB:
    def __init__(self, env, fob_id, convoys):
        self.env = env
        self.id = fob_id
        self.inventory = {supply: 100 for supply in SUPPLY_TYPES}  # Start at 100%
        self.convoys = convoys
        self.critical_level = 30  # Below 30% is critical
        
    def consume_supplies(self):
        """Consume supplies over time."""
        while True:
            yield self.env.timeout(24)  # Daily consumption
            
            for supply in SUPPLY_TYPES:
                # Variable consumption based on operational tempo
                consumption = DAILY_CONSUMPTION_RATE * random.uniform(0.8, 1.2)
                self.inventory[supply] = max(0, self.inventory[supply] - consumption * 100)
                
                # Check if resupply needed
                if self.inventory[supply] < self.critical_level:
                    self.env.process(self.request_resupply(supply))
            
            # Record status
            supply_data.append({
                'time': self.env.now,
                'fob_id': self.id,
                'ammunition': self.inventory['Ammunition'],
                'fuel': self.inventory['Fuel'],
                'food': self.inventory['Food'],
                'medical': self.inventory['Medical']
            })
    
    def request_resupply(self, supply_type):
        """Request emergency resupply for critical items."""
        priority = 1 if self.inventory[supply_type] < 15 else 2
        yield self.env.process(convoy_mission(self.env, self.id, supply_type, 
                                            self.convoys, priority))

def convoy_mission(env, fob_id, supply_type, convoys, priority):
    """Execute a supply convoy mission."""
    mission_start = env.now
    
    with convoys.request(priority=priority) as request:
        yield request
        convoy_wait = env.now - mission_start
        
        # Travel to FOB and deliver
        travel_time = random.gauss(DELIVERY_TIME_MEAN, DELIVERY_TIME_MEAN * 0.2)
        travel_time = max(12, travel_time)  # Minimum 12 hours
        yield env.timeout(travel_time)
        
        # Deliver supplies
        delivery_amount = random.uniform(40, 60)  # Resupply 40-60%
        
        # Update FOB inventory
        for fob in fobs:
            if fob.id == fob_id:
                old_level = fob.inventory[supply_type]
                fob.inventory[supply_type] = min(100, old_level + delivery_amount)
                
                convoy_missions.append({
                    'mission_time': mission_start,
                    'fob_id': fob_id,
                    'supply_type': supply_type,
                    'priority': priority,
                    'wait_time': convoy_wait,
                    'travel_time': travel_time,
                    'delivered': delivery_amount,
                    'old_level': old_level,
                    'new_level': fob.inventory[supply_type]
                })
                
                priority_name = 'URGENT' if priority == 1 else 'Regular'
                print(f"[{env.now/24:.1f}d] {priority_name} convoy delivered {delivery_amount:.0f}% "
                      f"{supply_type} to FOB-{fob_id} ({old_level:.0f}% → {fob.inventory[supply_type]:.0f}%)")
                break

# Run simulation
env = simpy.Environment()
convoys = simpy.PriorityResource(env, capacity=NUM_CONVOYS)
fobs = [FOB(env, i, convoys) for i in range(NUM_FOBS)]

# Start consumption processes
for fob in fobs:
    env.process(fob.consume_supplies())

env.run(until=SIM_TIME)

# Analysis
df_supply = pd.DataFrame(supply_data)
df_missions = pd.DataFrame(convoy_missions)

print(f"\n=== Supply Chain Performance ===")
print(f"Total convoy missions: {len(df_missions)}")
print(f"Emergency missions: {len(df_missions[df_missions['priority'] == 1])}")
print(f"Average response time: {df_missions['wait_time'].mean():.1f} hours")
print(f"Average delivery time: {df_missions['travel_time'].mean():.1f} hours")

# Calculate stockout incidents
stockouts = 0
for supply in SUPPLY_TYPES:
    stockouts += (df_supply[supply.lower()] == 0).sum()
print(f"Total stockout incidents: {stockouts}")

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Inventory levels over time for all FOBs
for fob_id in range(NUM_FOBS):
    fob_data = df_supply[df_supply['fob_id'] == fob_id]
    axes[0, 0].plot(fob_data['time']/24, 
                    (fob_data['ammunition'] + fob_data['fuel'] + 
                     fob_data['food'] + fob_data['medical'])/4,
                    label=f'FOB-{fob_id}', alpha=0.7, linewidth=2)
axes[0, 0].axhline(y=30, color='red', linestyle='--', alpha=0.5, label='Critical level')
axes[0, 0].set_xlabel('Time (days)')
axes[0, 0].set_ylabel('Average Inventory Level (%)')
axes[0, 0].set_title('FOB Inventory Levels Over Time')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. Supply type breakdown for last time period
last_time = df_supply['time'].max()
last_data = df_supply[df_supply['time'] == last_time]
supply_means = last_data[['ammunition', 'fuel', 'food', 'medical']].mean()
axes[0, 1].bar(SUPPLY_TYPES, supply_means.values, alpha=0.7, 
               color=['red', 'orange', 'green', 'blue'])
axes[0, 1].axhline(y=30, color='red', linestyle='--', alpha=0.5)
axes[0, 1].set_ylabel('Inventory Level (%)')
axes[0, 1].set_title('Current Inventory by Supply Type')
axes[0, 1].set_ylim(0, 100)
for i, v in enumerate(supply_means.values):
    axes[0, 1].text(i, v + 2, f'{v:.0f}%', ha='center')

# 3. Convoy utilization over time
mission_times = df_missions['mission_time'] / 24  # Convert to days
axes[1, 0].hist(mission_times, bins=30, alpha=0.7, color='purple', edgecolor='black')
axes[1, 0].set_xlabel('Time (days)')
axes[1, 0].set_ylabel('Number of Missions')
axes[1, 0].set_title('Convoy Mission Distribution Over Time')
axes[1, 0].grid(True, alpha=0.3, axis='y')

# 4. Mission priority and response times
priority_data = df_missions.groupby('priority').agg({
    'wait_time': 'mean',
    'travel_time': 'mean'
}).reset_index()
priority_data['priority_name'] = priority_data['priority'].map({1: 'Emergency', 2: 'Regular'})

x = np.arange(len(priority_data))
width = 0.35
axes[1, 1].bar(x - width/2, priority_data['wait_time'], width, 
               label='Wait Time', alpha=0.7, color='coral')
axes[1, 1].bar(x + width/2, priority_data['travel_time'], width, 
               label='Travel Time', alpha=0.7, color='skyblue')
axes[1, 1].set_xlabel('Mission Priority')
axes[1, 1].set_ylabel('Time (hours)')
axes[1, 1].set_title('Response Times by Mission Priority')
axes[1, 1].set_xticks(x)
axes[1, 1].set_xticklabels(priority_data['priority_name'])
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")  # Special marker for plot data

Output

Click "Run Simulation" to see results...

Building Climate Control

This simulation models an HVAC system managing climate control for a multi-zone building. The simulation demonstrates:

Zone temperature regulation and control
Energy consumption optimization
Occupancy-based climate adjustment
Peak load management and efficiency

Adjust parameters to balance comfort levels with energy efficiency across different building zones.

Simulation Parameters

Number of Zones

Number of HVAC Units

Target Temperature (°C)

Temperature Tolerance (°C)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_ZONES = {NUM_ZONES}
NUM_HVAC_UNITS = {NUM_HVAC_UNITS}
TARGET_TEMP = {TARGET_TEMP}
TEMP_TOLERANCE = {TEMP_TOLERANCE}
SIM_TIME = 1440  # 24 hours in minutes

# Zone characteristics
zone_data = []
hvac_operations = []
energy_consumption = []

class Zone:
    def __init__(self, env, zone_id, hvac_units):
        self.env = env
        self.id = zone_id
        self.hvac_units = hvac_units
        self.temperature = TARGET_TEMP + random.uniform(-3, 3)
        self.occupancy = 0
        self.heat_load = random.uniform(0.8, 1.2)  # Zone-specific heat characteristics
        self.size = random.uniform(80, 150)  # Square meters
        
    def temperature_drift(self):
        """Simulate temperature changes due to external factors."""
        while True:
            yield self.env.timeout(5)  # Check every 5 minutes
            
            # External temperature influence (simplified day/night cycle)
            hour = (self.env.now // 60) % 24
            external_temp = 15 + 10 * np.sin((hour - 6) * np.pi / 12)
            
            # Occupancy influence (more people = more heat)
            occupancy_heat = self.occupancy * 0.02
            
            # Temperature drift
            drift = (external_temp - self.temperature) * 0.05 * self.heat_load + occupancy_heat
            self.temperature += drift
            
            # Check if HVAC intervention needed
            if abs(self.temperature - TARGET_TEMP) > TEMP_TOLERANCE:
                self.env.process(self.request_hvac())
            
            # Record zone data
            zone_data.append({
                'time': self.env.now,
                'zone_id': self.id,
                'temperature': self.temperature,
                'occupancy': self.occupancy,
                'external_temp': external_temp,
                'size': self.size
            })
    
    def occupancy_pattern(self):
        """Simulate occupancy patterns throughout the day."""
        while True:
            hour = (self.env.now // 60) % 24
            
            # Office hours pattern (9 AM - 6 PM)
            if 9 <= hour <= 18:
                self.occupancy = random.randint(5, 20) * (self.size / 100)
            else:
                self.occupancy = random.randint(0, 3)
            
            yield self.env.timeout(60)  # Update hourly
    
    def request_hvac(self):
        """Request HVAC adjustment."""
        with self.hvac_units.request() as request:
            yield request
            
            # Calculate required adjustment
            temp_diff = self.temperature - TARGET_TEMP
            adjustment_time = abs(temp_diff) * 2  # Minutes to adjust
            energy_rate = abs(temp_diff) * self.size * 0.1  # kW
            
            operation_start = self.env.now
            
            # Perform temperature adjustment
            yield self.env.timeout(adjustment_time)
            self.temperature = TARGET_TEMP + random.uniform(-0.5, 0.5)
            
            # Record HVAC operation
            hvac_operations.append({
                'start_time': operation_start,
                'end_time': self.env.now,
                'zone_id': self.id,
                'temp_adjustment': temp_diff,
                'energy_used': energy_rate * adjustment_time / 60  # kWh
            })

def monitor_energy(env):
    """Monitor total building energy consumption."""
    while True:
        yield env.timeout(15)  # Every 15 minutes
        
        # Calculate energy usage in last 15 minutes
        current_time = env.now
        recent_ops = [op for op in hvac_operations 
                     if op['end_time'] >= current_time - 15]
        
        # Sum energy used in the last 15-minute window
        total_energy = sum(op['energy_used'] for op in recent_ops)
        # Convert to average power over 15 minutes (kW)
        current_power = total_energy * 4  # (kWh * 4 = kW for 15-min intervals)
        
        energy_consumption.append({
            'time': current_time,
            'power': current_power,
            'recent_operations': len(recent_ops)
        })

# Run simulation
env = simpy.Environment()
hvac_units = simpy.Resource(env, capacity=NUM_HVAC_UNITS)
zones = [Zone(env, i, hvac_units) for i in range(NUM_ZONES)]

# Start processes
for zone in zones:
    env.process(zone.temperature_drift())
    env.process(zone.occupancy_pattern())

env.process(monitor_energy(env))
env.run(until=SIM_TIME)

# Analysis
df_zones = pd.DataFrame(zone_data)
df_hvac = pd.DataFrame(hvac_operations)
df_energy = pd.DataFrame(energy_consumption)

print(f"\n=== HVAC System Performance ===")
print(f"Total HVAC operations: {len(df_hvac)}")
print(f"Total energy consumed: {df_hvac['energy_used'].sum():.1f} kWh")
print(f"Average zone temperature: {df_zones['temperature'].mean():.2f}°C")
print(f"Temperature variance: {df_zones['temperature'].std():.2f}°C")

# Comfort analysis
comfort_violations = ((df_zones['temperature'] < TARGET_TEMP - TEMP_TOLERANCE) | 
                     (df_zones['temperature'] > TARGET_TEMP + TEMP_TOLERANCE)).sum()
comfort_rate = (1 - comfort_violations / len(df_zones)) * 100
print(f"Comfort level maintained: {comfort_rate:.1f}% of time")

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Zone temperatures over time
for zone_id in range(min(4, NUM_ZONES)):  # Show first 4 zones
    zone_temps = df_zones[df_zones['zone_id'] == zone_id]
    axes[0, 0].plot(zone_temps['time']/60, zone_temps['temperature'], 
                    label=f'Zone {zone_id}', alpha=0.7)

axes[0, 0].axhline(y=TARGET_TEMP, color='green', linestyle='--', label='Target')
axes[0, 0].axhline(y=TARGET_TEMP + TEMP_TOLERANCE, color='red', linestyle=':', alpha=0.5)
axes[0, 0].axhline(y=TARGET_TEMP - TEMP_TOLERANCE, color='red', linestyle=':', alpha=0.5)
axes[0, 0].set_xlabel('Time (hours)')
axes[0, 0].set_ylabel('Temperature (°C)')
axes[0, 0].set_title('Zone Temperatures Throughout Day')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. Energy consumption profile
if len(df_energy) > 0:
    axes[0, 1].plot(df_energy['time']/60, df_energy['power'], 
                    linewidth=2, color='orange')
    axes[0, 1].fill_between(df_energy['time']/60, df_energy['power'], 
                           alpha=0.3, color='orange')
axes[0, 1].set_xlabel('Time (hours)')
axes[0, 1].set_ylabel('Power Consumption (kW)')
axes[0, 1].set_title('Building Energy Consumption Profile')
axes[0, 1].grid(True, alpha=0.3)

# 3. Temperature distribution by zone
zone_temp_avg = df_zones.groupby('zone_id')['temperature'].agg(['mean', 'std'])
axes[1, 0].bar(zone_temp_avg.index, zone_temp_avg['mean'], 
               yerr=zone_temp_avg['std'], alpha=0.7, capsize=5)
axes[1, 0].axhline(y=TARGET_TEMP, color='green', linestyle='--')
axes[1, 0].set_xlabel('Zone ID')
axes[1, 0].set_ylabel('Average Temperature (°C)')
axes[1, 0].set_title('Average Temperature by Zone')
axes[1, 0].grid(True, alpha=0.3, axis='y')

# 4. HVAC utilization and efficiency
hourly_ops = df_hvac.copy()
hourly_ops['hour'] = hourly_ops['start_time'] // 60
hourly_energy = hourly_ops.groupby('hour')['energy_used'].sum()
hourly_count = hourly_ops.groupby('hour').size()

ax4 = axes[1, 1]
ax4_twin = ax4.twinx()

bars = ax4.bar(hourly_energy.index, hourly_energy.values, alpha=0.7, 
                color='steelblue', label='Energy (kWh)')
line = ax4_twin.plot(hourly_count.index, hourly_count.values, 
                     color='red', marker='o', linewidth=2, label='Operations')

ax4.set_xlabel('Hour of Day')
ax4.set_ylabel('Energy Consumption (kWh)', color='steelblue')
ax4_twin.set_ylabel('Number of Operations', color='red')
ax4.set_title('HVAC Operations and Energy Use by Hour')
ax4.tick_params(axis='y', labelcolor='steelblue')
ax4_twin.tick_params(axis='y', labelcolor='red')
ax4.grid(True, alpha=0.3)

# Add legends
lines1, labels1 = ax4.get_legend_handles_labels()
lines2, labels2 = ax4_twin.get_legend_handles_labels()
ax4.legend(lines1 + lines2, labels1 + labels2, loc='upper left')

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")  # Special marker for plot data

Output

Click "Run Simulation" to see results...

Multi-Year Guaranteed Annuity (MYGA)

This simulation models a portfolio of Multi-Year Guaranteed Annuities (MYGAs). The simulation demonstrates:

Interest accumulation at a guaranteed rate
Policyholder behavior, including lapses (surrenders)
Cash flow analysis for an insurance company
Portfolio valuation over time

Adjust the parameters to see how interest rates and policyholder behavior impact the growth and stability of the annuity portfolio.

Simulation Parameters

Number of Policies

Initial Investment ($)

Annual Interest Rate (%)

Annuity Term (years)

Annual Lapse Rate (%)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_POLICIES = {NUM_POLICIES}
INITIAL_INVESTMENT_MEAN = {INITIAL_INVESTMENT_MEAN}
INTEREST_RATE = {INTEREST_RATE} / 100
ANNUITY_TERM = {ANNUITY_TERM}
LAPSE_RATE = {LAPSE_RATE} / 100
SIM_TIME = ANNUITY_TERM + 1

# Data collection
policy_data = []
portfolio_value_over_time = []

def annuity_policy(env, policy_id):
    """Simulate a single annuity policy."""
    initial_investment = random.uniform(INITIAL_INVESTMENT_MEAN * 0.8, INITIAL_INVESTMENT_MEAN * 1.2)
    account_value = initial_investment
    
    policy_data.append({
        'policy_id': policy_id,
        'year': 0,
        'start_value': account_value,
        'interest': 0,
        'end_value': account_value,
        'status': 'Active'
    })
    
    for year in range(1, ANNUITY_TERM + 1):
        yield env.timeout(1) # Simulate one year
        
        # Check for lapse
        if random.random() < LAPSE_RATE:
            policy_data.append({
                'policy_id': policy_id,
                'year': year,
                'start_value': account_value,
                'interest': 0,
                'end_value': account_value,
                'status': 'Lapsed'
            })
            break
            
        # Accrue interest
        interest_earned = account_value * INTEREST_RATE
        account_value += interest_earned
        
        policy_data.append({
            'policy_id': policy_id,
            'year': year,
            'start_value': account_value - interest_earned,
            'interest': interest_earned,
            'end_value': account_value,
            'status': 'Active'
        })
    else: # If loop completed without break
        # Update status to Matured
        policy_data[-1]['status'] = 'Matured'

def policy_generator(env):
    """Generate all policies at the start of the simulation."""
    for i in range(NUM_POLICIES):
        env.process(annuity_policy(env, i))

# Run simulation
env = simpy.Environment()
policy_generator(env)
env.run(until=SIM_TIME)

# Analysis
df = pd.DataFrame(policy_data)

print(f"\n=== Annuity Portfolio Analysis ===")
final_year_df = df[df['year'] == ANNUITY_TERM]
matured_policies = final_year_df[final_year_df['status'] == 'Matured']
lapsed_policies = df[df['status'] == 'Lapsed']

print(f"Total policies simulated: {NUM_POLICIES}")
print(f"Matured policies: {len(matured_policies)}")
print(f"Lapsed policies: {len(lapsed_policies)}")

if not matured_policies.empty:
    avg_final_value = matured_policies['end_value'].mean()
    print(f"Average matured account value: ${avg_final_value:,.2f}")
    total_interest = df[df['policy_id'].isin(matured_policies['policy_id'])]['interest'].sum()
    print(f"Total interest paid on matured policies: ${total_interest:,.2f}")

if not lapsed_policies.empty:
    avg_lapse_year = lapsed_policies['year'].mean()
    print(f"Average policy year at lapse: {avg_lapse_year:.1f}")

# Portfolio value over time
portfolio_value = df.groupby('year')['end_value'].sum()

# Create visualizations
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Portfolio Value Over Time
axes[0, 0].plot(portfolio_value.index, portfolio_value.values / 1e6, marker='o', linewidth=2)
axes[0, 0].set_xlabel('Year')
axes[0, 0].set_ylabel('Portfolio Value ($ Millions)')
axes[0, 0].set_title('Total Portfolio Value Over Time')
axes[0, 0].grid(True, alpha=0.3)

# 2. Distribution of Final Account Values (Matured Policies)
if not matured_policies.empty:
    sns.histplot(matured_policies['end_value'], bins=30, kde=True, ax=axes[0, 1])
    axes[0, 1].set_xlabel('Final Account Value ($)')
    axes[0, 1].set_ylabel('Number of Policies')
    axes[0, 1].set_title('Distribution of Matured Policy Values')

# 3. Active vs. Lapsed Policies Over Time
active_counts = df[df['status'] != 'Lapsed'].groupby('year')['policy_id'].nunique()
lapsed_counts = df[df['status'] == 'Lapsed'].groupby('year')['policy_id'].nunique().cumsum()

axes[1, 0].stackplot(active_counts.index, active_counts.values, lapsed_counts.reindex(active_counts.index, fill_value=0).values,
                     labels=['Active', 'Lapsed'],
                     colors=['#3498db', '#e74c3c'], alpha=0.7)
axes[1, 0].set_xlabel('Year')
axes[1, 0].set_ylabel('Number of Policies')
axes[1, 0].set_title('Policy Status Over Time')
axes[1, 0].legend(loc='upper left')
axes[1, 0].grid(True, alpha=0.3)

# 4. Total Interest Earned Per Year
yearly_interest = df.groupby('year')['interest'].sum()
axes[1, 1].bar(yearly_interest.index, yearly_interest.values / 1e6, alpha=0.7)
axes[1, 1].set_xlabel('Year')
axes[1, 1].set_ylabel('Total Interest ($ Millions)')
axes[1, 1].set_title('Total Interest Earned by Year')
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")

Output

Click "Run Simulation" to see results...

Actuarial Pricing

This simulation models the process an actuary uses to price an insurance premium for a portfolio of policies. It uses a Monte Carlo approach to simulate one year of claims activity many times to understand the range of possible outcomes. The simulation demonstrates:

Modeling claim frequency and severity
Calculating pure premium based on expected losses
Applying a loading factor for expenses and profit
Analyzing the risk of the insurance portfolio

Adjust the parameters to see how claim characteristics and business assumptions affect the final premium and the insurer's risk.

Simulation Parameters

Number of Policies

Annual Claim Frequency (per policy)

Mean Claim Severity ($)

Claim Severity Std Dev ($)

Expense & Profit Loading (%)

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_POLICIES = {NUM_POLICIES}
CLAIM_FREQUENCY = {CLAIM_FREQUENCY}  # per policy per year
SEVERITY_MEAN = {SEVERITY_MEAN}
SEVERITY_STD = {SEVERITY_STD}
LOADING_FACTOR = {LOADING_FACTOR} / 100
NUM_SIMULATIONS = 100  # Monte Carlo runs
SIM_DURATION = 1  # Simulate for 1 year

# Data collection
total_claims_per_run = []

def get_lognorm_params(mean, stddev):
    """Convert mean and stddev to lognormal distribution parameters."""
    if mean <= 0: return -1, 1
    variance = stddev**2
    mu = np.log(mean**2 / np.sqrt(variance + mean**2))
    sigma = np.sqrt(np.log(variance/mean**2 + 1))
    return mu, sigma

# Calculate lognormal parameters for claim severity
SEVERITY_MU, SEVERITY_SIGMA = get_lognorm_params(SEVERITY_MEAN, SEVERITY_STD)

def claim_generator(env, claims_list):
    """Generate claims over the simulation period."""
    total_expected_claims = NUM_POLICIES * CLAIM_FREQUENCY
    if total_expected_claims == 0:
        return
        
    while True:
        # Time to next claim follows an exponential distribution
        yield env.timeout(random.expovariate(total_expected_claims))
        
        # Generate claim severity from a log-normal distribution
        claim_amount = random.lognormvariate(SEVERITY_MU, SEVERITY_SIGMA)
        claims_list.append(claim_amount)

# Run Monte Carlo simulation
for i in range(NUM_SIMULATIONS):
    env = simpy.Environment()
    claims_this_run = []
    env.process(claim_generator(env, claims_this_run))
    env.run(until=SIM_DURATION)
    
    total_claims = sum(claims_this_run)
    total_claims_per_run.append(total_claims)

# Analysis
df_results = pd.DataFrame({'total_claims': total_claims_per_run})

# Premium Calculation
expected_total_claims = df_results['total_claims'].mean()
pure_premium_per_policy = expected_total_claims / NUM_POLICIES
gross_premium_per_policy = pure_premium_per_policy * (1 + LOADING_FACTOR)
total_premium_income = gross_premium_per_policy * NUM_POLICIES

# Risk Metrics
loss_ratio = expected_total_claims / total_premium_income * 100
probability_of_loss = (df_results['total_claims'] > total_premium_income).mean() * 100

print(f"=== Actuarial Pricing Simulation Results ({NUM_SIMULATIONS} runs) ===")
print(f"\n--- Premium Calculation ---")
print(f"Expected Annual Claims (Portfolio): ${expected_total_claims:,.2f}")
print(f"Pure Premium per Policy: ${pure_premium_per_policy:,.2f}")
print(f"Gross Premium per Policy (w/ loading): ${gross_premium_per_policy:,.2f}")
print(f"Total Premium Income (Portfolio): ${total_premium_income:,.2f}")

print(f"\n--- Risk Analysis ---")
print(f"Projected Loss Ratio: {loss_ratio:.1f}%")
print(f"Probability of Annual Loss (Claims > Premiums): {probability_of_loss:.1f}%")
print(f"Standard Deviation of Annual Claims: ${df_results['total_claims'].std():,.2f}")
print(f"95th Percentile of Annual Claims: ${df_results['total_claims'].quantile(0.95):,.2f}")

# Create visualizations
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 1. Distribution of Total Annual Claims
sns.histplot(df_results['total_claims'], bins=30, kde=True, ax=axes[0])
axes[0].axvline(expected_total_claims, color='green', linestyle='--', label=f'Mean: ${expected_total_claims:,.0f}')
axes[0].axvline(total_premium_income, color='red', linestyle='-', label=f'Premium Income: ${total_premium_income:,.0f}')
axes[0].set_xlabel('Total Annual Claims ($)')
axes[0].set_ylabel('Frequency')
axes[0].set_title(f'Distribution of Total Claims ({NUM_SIMULATIONS} runs)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 2. Example Claim Severity Distribution
example_claims = [random.lognormvariate(SEVERITY_MU, SEVERITY_SIGMA) for _ in range(1000)]
sns.histplot(example_claims, bins=30, kde=True, ax=axes[1], color='purple')
axes[1].axvline(np.mean(example_claims), color='orange', linestyle='--', label=f'Sample Mean: ${np.mean(example_claims):,.0f}')
axes[1].set_xlabel('Individual Claim Severity ($)')
axes[1].set_ylabel('Frequency')
axes[1].set_title('Example Individual Claim Severity')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")

Output

Click "Run Simulation" to see results...

Actuarial Reserving

This simulation models the process of establishing reserves for insurance claims. It simulates claims occurring in a single "accident year" and tracks how they are reported and paid over subsequent "development years". This helps estimate IBNR (Incurred But Not Reported) reserves. The simulation demonstrates:

Modeling reporting and payment delays for claims
Tracking the development of paid losses over time
Estimating ultimate losses for an accident year
Calculating the required reserves (Outstanding Claims)

Adjust the parameters to see how claim delays and other factors influence the loss development pattern and the size of the required reserves.

Simulation Parameters

Number of Policies

Annual Claim Frequency

Mean Claim Severity ($)

Mean Reporting Lag (Years)

Mean Payment Lag (Years)

Development Years

Simulation Code

import simpy
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import io

# Set random seed for reproducibility
random.seed(42)
np.random.seed(42)

# Simulation parameters
NUM_POLICIES = {NUM_POLICIES}
CLAIM_FREQUENCY = {CLAIM_FREQUENCY}
SEVERITY_MEAN = {SEVERITY_MEAN}
REPORTING_LAG_MEAN = {REPORTING_LAG_MEAN}
PAYMENT_LAG_MEAN = {PAYMENT_LAG_MEAN}  # After reporting
DEVELOPMENT_YEARS = {DEVELOPMENT_YEARS}

# Data collection
claims_data = []

def claim_lifecycle(env, claim_id, severity):
    """Simulates a single claim's journey from occurrence to payment."""
    occurrence_time = env.now
    
    # Simulate reporting and payment delays
    reporting_delay = random.expovariate(1.0 / REPORTING_LAG_MEAN)
    payment_delay = random.expovariate(1.0 / PAYMENT_LAG_MEAN)
    
    # Wait for the claim to be reported
    yield env.timeout(reporting_delay)
    report_time = env.now
    
    # Wait for the claim to be paid
    yield env.timeout(payment_delay)
    payment_time = env.now
    
    claims_data.append({
        'claim_id': claim_id,
        'severity': severity,
        'occurrence_time': occurrence_time,
        'report_time': report_time,
        'payment_time': payment_time
    })

def accident_year_generator(env):
    """Generates all claims for a single accident year (Year 0)."""
    num_claims = np.random.poisson(NUM_POLICIES * CLAIM_FREQUENCY)
    
    # To model claims occurring uniformly, we generate their occurrence times and sort them
    occurrence_times = sorted([random.uniform(0, 1) for _ in range(num_claims)])

for i, occurrence_time in enumerate(occurrence_times):
        # Wait until the next scheduled occurrence time. This prevents negative timeouts.
        yield env.timeout(occurrence_time - env.now)
        
        severity = random.expovariate(1.0 / SEVERITY_MEAN)
        env.process(claim_lifecycle(env, i, severity))

# Run simulation
env = simpy.Environment()
env.process(accident_year_generator(env))
env.run(until=DEVELOPMENT_YEARS)

# Analysis
df = pd.DataFrame(claims_data)
ultimate_loss = df['severity'].sum()

# Calculate paid amounts at the end of each development year
paid_by_year = []
for year in range(DEVELOPMENT_YEARS + 1):
    paid = df[df['payment_time'] <= year]['severity'].sum()
    paid_by_year.append({'year': year, 'paid': paid})

df_paid = pd.DataFrame(paid_by_year)
df_paid['reserve'] = ultimate_loss - df_paid['paid']

print("=== Actuarial Reserving Simulation Results ===")
print(f"\nAccident Year Summary:")
print(f"Total Claims Incurred: {len(df)}")
print(f"True Ultimate Loss: ${ultimate_loss:,.2f}")

print("\nLoss Development:")
for _, row in df_paid.iterrows():
    if row['year'] <= DEVELOPMENT_YEARS:
        print(f"Year {row['year']}: Paid=${row['paid']:,.0f}, Reserve=${row['reserve']:,.0f}")

# Create visualizations
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 1. Cumulative Paid Claims (Loss Development)
axes[0].plot(df_paid['year'], df_paid['paid'], marker='o', label='Cumulative Paid')
axes[0].axhline(y=ultimate_loss, color='red', linestyle='--', label='Ultimate Loss')
axes[0].set_xlabel('Development Year')
axes[0].set_ylabel('Amount ($)')
axes[0].set_title('Loss Development Over Time')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(bottom=0)

# 2. Reserve vs. Paid at each year
width = 0.8
axes[1].bar(df_paid['year'], df_paid['paid'], width, label='Paid to Date', color='skyblue')
axes[1].bar(df_paid['year'], df_paid['reserve'], width, bottom=df_paid['paid'], 
            label='Outstanding Reserve', color='coral')
axes[1].set_xlabel('Development Year')
axes[1].set_ylabel('Amount ($)')
axes[1].set_title('Paid vs. Reserve by Year')
axes[1].axhline(y=ultimate_loss, color='red', linestyle='--', alpha=0.7)
axes[1].legend()
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()

# Create a base64 encoded plot for web display
import base64
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
plot_base64 = base64.b64encode(buf.read()).decode()
buf.close()
print(f"PLOT_DATA:{plot_base64}")

Output

Click "Run Simulation" to see results...

SimPy Simulation Playground

Steel Production Plant Simulation

Simulation Parameters

Simulation Code

Output

Shipping Logistics Simulation

Simulation Parameters

Simulation Code

Output

Satellite Communication Network

Simulation Parameters

Simulation Code

Output

Aircraft Turnaround Operations

Simulation Parameters

Simulation Code

Output

Military Supply Chain

Simulation Parameters

Simulation Code

Output

Building Climate Control

Simulation Parameters

Simulation Code

Output

Multi-Year Guaranteed Annuity (MYGA)

Simulation Parameters

Simulation Code

Output

Actuarial Pricing

Simulation Parameters

Simulation Code

Output

Actuarial Reserving

Simulation Parameters

Simulation Code

Output

Would you like a copy of my book on Simulation in Python with SimPy?