Navigating Uncertainty：Solving the Stochastic Vehicle Routing Problem with Python

May 7, 2025

In today’s fast-paced logistics world, planning efficient delivery routes isn’t just about finding the shortest path between points.
Real-world delivery operations face constant uncertainty - customer demands fluctuate, traffic conditions change, and unexpected events disrupt even the best-laid plans.

This is where the Stochastic Vehicle Routing Problem (SVRP) comes into play.
Unlike its deterministic counterpart, SVRP acknowledges that customer demands aren’t fixed but follow probability distributions.
Today, I’ll walk through a practical example of SVRP and solve it using Python.

Understanding the Stochastic Vehicle Routing Problem

In the SVRP, we need to design optimal routes for a fleet of vehicles to serve customers whose demands are uncertain.
The objective is typically to minimize the expected total cost while ensuring that the probability of route failure (where demand exceeds vehicle capacity) stays below an acceptable threshold.

The mathematical formulation of SVRP can be expressed as:

$$\min E[f(x,\xi)]$$

where $x$ represents the routing decisions, $\xi$ represents the random demand, and $f(x,\xi)$ is the cost function.
The constraints include:

$$P(g_i(x,\xi) \leq 0) \geq 1-\alpha_i, \forall i$$

Here, $g_i$ represents constraints (like capacity constraints), and $\alpha_i$ is the acceptable risk level.

Building a Practical SVRP Example

Let’s tackle a concrete example: A company operates a distribution center that needs to deliver products to 10 customers.
The demands of these customers are not fixed but follow normal distributions.
The company wants to determine optimal routes that minimize the total distance traveled while keeping the risk of exceeding vehicle capacity below 10%.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from scipy.stats import norm
from scipy.spatial import distance_matrix
import random
from itertools import permutations
import networkx as nx
from matplotlib.patches import Patch

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

# Problem parameters
num_customers = 10
num_vehicles = 3
vehicle_capacity = 30
depot_coords = (0, 0)
confidence_level = 0.9  # 90% confidence that routes won't fail due to capacity

# Generate customer locations
customer_x = np.random.uniform(-10, 10, num_customers)
customer_y = np.random.uniform(-10, 10, num_customers)
customer_coords = list(zip(customer_x, customer_y))

# Add depot as node 0
all_coords = [depot_coords] + customer_coords

# Generate stochastic demands (normal distribution)
demand_means = np.random.uniform(5, 15, num_customers)
demand_stds = demand_means * 0.3  # Standard deviation as 30% of mean

# Calculate distance matrix
dist_matrix = np.zeros((num_customers + 1, num_customers + 1))
for i in range(num_customers + 1):
    for j in range(num_customers + 1):
        if i != j:
            dist_matrix[i, j] = np.sqrt((all_coords[i][0] - all_coords[j][0])**2 +
                                        (all_coords[i][1] - all_coords[j][1])**2)

# Function to calculate route length
def calculate_route_length(route, dist_mat):
    length = 0
    prev_node = 0  # Start at depot
    for node in route:
        length += dist_mat[prev_node, node]
        prev_node = node
    length += dist_mat[prev_node, 0]  # Return to depot
    return length

# Function to check if a route satisfies capacity constraints with given confidence
def check_capacity_constraint(route, means, stds, capacity, confidence):
    route_mean_demand = sum(means[i-1] for i in route)
    route_var_demand = sum(stds[i-1]**2 for i in route)
    route_std_demand = np.sqrt(route_var_demand)

    # Calculate the critical value based on confidence level
    z_value = norm.ppf(confidence)

    # The route will satisfy capacity with the given confidence if:
    # mean + z * std <= capacity
    effective_demand = route_mean_demand + z_value * route_std_demand

    return effective_demand <= capacity, effective_demand

# Savings algorithm for SVRP
def savings_algorithm_svrp(dist_mat, means, stds, capacity, confidence, num_nodes):
    # Calculate savings for each pair of customers
    savings = {}
    for i in range(1, num_nodes):
        for j in range(1, num_nodes):
            if i != j:
                # Savings formula: s(i,j) = d(0,i) + d(0,j) - d(i,j)
                savings[(i, j)] = dist_mat[0, i] + dist_mat[0, j] - dist_mat[i, j]

    # Sort savings in descending order
    sorted_savings = sorted(savings.items(), key=lambda x: x[1], reverse=True)

    # Initialize routes: one route per customer
    routes = [[i] for i in range(1, num_nodes)]

    # Try to merge routes based on savings
    for (i, j), saving in sorted_savings:
        # Find routes containing i and j
        route_i = None
        route_j = None
        for r in routes:
            if i in r and i == r[-1]:
                route_i = r
            if j in r and j == r[0]:
                route_j = r

        # If i and j are in different routes and at the right positions
        if route_i and route_j and route_i != route_j:
            # Check if merging would violate capacity constraints
            merged_route = route_i + route_j
            feasible, _ = check_capacity_constraint(merged_route, means, stds, capacity, confidence)

            if feasible:
                # Merge routes
                routes.remove(route_i)
                routes.remove(route_j)
                routes.append(merged_route)

    return routes

# Apply the savings algorithm
routes = savings_algorithm_svrp(dist_matrix, demand_means, demand_stds,
                               vehicle_capacity, confidence_level, num_customers + 1)

# Evaluate the solution
total_distance = 0
route_demands = []
route_risks = []

for route in routes:
    route_length = calculate_route_length(route, dist_matrix)
    total_distance += route_length

    # Calculate route statistics
    feasible, effective_demand = check_capacity_constraint(route, demand_means, demand_stds,
                                                          vehicle_capacity, confidence_level)
    route_mean_demand = sum(demand_means[i-1] for i in route)
    route_std_demand = np.sqrt(sum(demand_stds[i-1]**2 for i in route))

    # Calculate risk of exceeding capacity
    risk = 1 - norm.cdf((vehicle_capacity - route_mean_demand) / route_std_demand)

    route_demands.append((route_mean_demand, route_std_demand, effective_demand))
    route_risks.append(risk)

# Output results
print(f"Number of routes: {len(routes)}")
print(f"Total expected distance: {total_distance:.2f}")

# Create dataframe for route statistics
route_stats = []
for i, route in enumerate(routes):
    mean_demand, std_demand, effective_demand = route_demands[i]
    route_stats.append({
        'Route': i+1,
        'Customers': route,
        'Mean Demand': mean_demand,
        'Std Demand': std_demand,
        'Effective Demand': effective_demand,
        'Risk of Exceeding Capacity': route_risks[i] * 100,  # Convert to percentage
        'Distance': calculate_route_length(route, dist_matrix)
    })

route_stats_df = pd.DataFrame(route_stats)
print("\nRoute Statistics:")
print(route_stats_df[['Route', 'Customers', 'Mean Demand', 'Effective Demand',
                      'Risk of Exceeding Capacity', 'Distance']])

# Visualize the solution
plt.figure(figsize=(12, 10))

# Plot customers and depot
plt.scatter(depot_coords[0], depot_coords[1], c='red', s=200, marker='*', label='Depot')
plt.scatter(customer_x, customer_y, c='blue', s=100, label='Customers')

# Add customer labels
for i, (x, y) in enumerate(customer_coords):
    plt.annotate(f"{i+1}", (x, y), fontsize=12)

# Define colors for routes
colors = ['green', 'purple', 'orange', 'brown', 'pink', 'gray', 'olive', 'cyan']

# Plot routes
for i, route in enumerate(routes):
    color = colors[i % len(colors)]

    # Plot route from depot to first customer
    plt.plot([depot_coords[0], customer_coords[route[0]-1][0]],
             [depot_coords[1], customer_coords[route[0]-1][1]],
             c=color, linewidth=2)

    # Plot route between customers
    for j in range(len(route)-1):
        plt.plot([customer_coords[route[j]-1][0], customer_coords[route[j+1]-1][0]],
                 [customer_coords[route[j]-1][1], customer_coords[route[j+1]-1][1]],
                 c=color, linewidth=2)

    # Plot route from last customer back to depot
    plt.plot([customer_coords[route[-1]-1][0], depot_coords[0]],
             [customer_coords[route[-1]-1][1], depot_coords[1]],
             c=color, linewidth=2,
             label=f"Route {i+1}: {route}")

plt.title('Stochastic Vehicle Routing Problem Solution', fontsize=16)
plt.xlabel('X Coordinate', fontsize=14)
plt.ylabel('Y Coordinate', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(loc='best', fontsize=12)
plt.axis('equal')

# Generate legend for routes
route_patches = []
for i, route in enumerate(routes):
    color = colors[i % len(colors)]
    route_patches.append(Patch(color=color, label=f"Route {i+1}: {route}"))

plt.legend(handles=[
    plt.Line2D([0], [0], marker='*', color='w', markerfacecolor='red', markersize=15, label='Depot'),
    plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='blue', markersize=10, label='Customers')
] + route_patches, loc='best', fontsize=12)

plt.tight_layout()
plt.savefig('svrp_solution.png', dpi=300)
plt.show()

# Visualize demand uncertainty
plt.figure(figsize=(14, 8))

# Plot demand distributions for each customer
x = np.linspace(0, 25, 100)
for i in range(num_customers):
    y = norm.pdf(x, demand_means[i], demand_stds[i])
    plt.plot(x, y, label=f"Customer {i+1}")

plt.axvline(x=vehicle_capacity, color='r', linestyle='--', label='Vehicle Capacity')
plt.title('Customer Demand Distributions', fontsize=16)
plt.xlabel('Demand', fontsize=14)
plt.ylabel('Probability Density', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('demand_distributions.png', dpi=300)
plt.show()

# Visualize route demand uncertainty
plt.figure(figsize=(12, 8))

# Plot cumulative route demand distributions
x = np.linspace(0, 60, 100)
for i, route in enumerate(routes):
    mean_demand, std_demand, _ = route_demands[i]
    y = norm.pdf(x, mean_demand, std_demand)
    plt.plot(x, y, label=f"Route {i+1}: {route}")

plt.axvline(x=vehicle_capacity, color='r', linestyle='--', label='Vehicle Capacity')
plt.title('Route Demand Distributions', fontsize=16)
plt.xlabel('Total Route Demand', fontsize=14)
plt.ylabel('Probability Density', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('route_demand_distributions.png', dpi=300)
plt.show()

# Risk analysis visualization
plt.figure(figsize=(10, 6))
risk_df = pd.DataFrame({
    'Route': [f"Route {i+1}" for i in range(len(routes))],
    'Risk (%)': [r * 100 for r in route_risks]
})
sns.barplot(x='Route', y='Risk (%)', data=risk_df)
plt.axhline(y=(1-confidence_level)*100, color='r', linestyle='--',
            label=f'Acceptable Risk ({(1-confidence_level)*100}%)')
plt.title('Risk of Exceeding Vehicle Capacity by Route', fontsize=16)
plt.xlabel('Route', fontsize=14)
plt.ylabel('Risk (%)', fontsize=14)
plt.grid(True, axis='y', linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('route_risks.png', dpi=300)
plt.show()

# Monte Carlo simulation to validate solution
def simulate_demands(means, stds, num_simulations=1000):
    return np.random.normal(
        np.tile(means, (num_simulations, 1)),
        np.tile(stds, (num_simulations, 1))
    )

# Number of simulations
num_simulations = 10000

# Simulate demands for all customers
simulated_demands = simulate_demands(demand_means, demand_stds, num_simulations)

# Calculate route failures
route_failures = np.zeros(len(routes))
for sim in range(num_simulations):
    for i, route in enumerate(routes):
        # Sum the demands for this route in this simulation
        route_demand = sum(max(0, simulated_demands[sim, j-1]) for j in route)
        if route_demand > vehicle_capacity:
            route_failures[i] += 1

# Calculate empirical failure probabilities
empirical_risks = route_failures / num_simulations * 100

# Compare theoretical and empirical risks
risk_comparison = pd.DataFrame({
    'Route': [f"Route {i+1}" for i in range(len(routes))],
    'Theoretical Risk (%)': [r * 100 for r in route_risks],
    'Empirical Risk (%)': empirical_risks
})

plt.figure(figsize=(12, 6))
risk_comparison_melted = pd.melt(risk_comparison, id_vars=['Route'],
                                var_name='Risk Type', value_name='Risk (%)')
sns.barplot(x='Route', y='Risk (%)', hue='Risk Type', data=risk_comparison_melted)
plt.axhline(y=(1-confidence_level)*100, color='r', linestyle='--',
            label=f'Acceptable Risk ({(1-confidence_level)*100}%)')
plt.title('Theoretical vs. Empirical Risk Comparison', fontsize=16)
plt.xlabel('Route', fontsize=14)
plt.ylabel('Risk (%)', fontsize=14)
plt.grid(True, axis='y', linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('risk_comparison.png', dpi=300)
plt.show()

print("\nRisk Comparison (Theoretical vs. Empirical):")
print(risk_comparison)

Code Implementation Explained

Let’s break down our implementation by key components:

1. Problem Setup

First, we set up our problem parameters and generate random customer locations and stochastic demands:

# Problem parameters
num_customers = 10
num_vehicles = 3
vehicle_capacity = 30
depot_coords = (0, 0)
confidence_level = 0.9  # 90% confidence that routes won't fail due to capacity

# Generate customer locations
customer_x = np.random.uniform(-10, 10, num_customers)
customer_y = np.random.uniform(-10, 10, num_customers)
customer_coords = list(zip(customer_x, customer_y))

# Generate stochastic demands (normal distribution)
demand_means = np.random.uniform(5, 15, num_customers)
demand_stds = demand_means * 0.3  # Standard deviation as 30% of mean

Here, we’re creating 10 customers with randomly generated coordinates.
The demands follow normal distributions with means between 5 and 15 units, and standard deviations at 30% of the respective means.
This reflects real-world scenarios where customer demands fluctuate around expected values.

2. Handling Demand Uncertainty

The core of SVRP is managing uncertainty.
For this, we implement a function that checks if routes satisfy capacity constraints with a given confidence level:

def check_capacity_constraint(route, means, stds, capacity, confidence):
    route_mean_demand = sum(means[i-1] for i in route)
    route_var_demand = sum(stds[i-1]**2 for i in route)
    route_std_demand = np.sqrt(route_var_demand)
    
    # Calculate the critical value based on confidence level
    z_value = norm.ppf(confidence)
    
    # The route will satisfy capacity with the given confidence if:
    # mean + z * std <= capacity
    effective_demand = route_mean_demand + z_value * route_std_demand
    
    return effective_demand <= capacity, effective_demand

This function applies probability theory to calculate an “effective demand” that represents the demand level we’re confident won’t be exceeded (in this case, with 90% confidence).
It uses the properties of normal distributions where:

$$P(X \leq \mu + z \cdot \sigma) = \Phi(z)$$

Where $\Phi$ is the cumulative distribution function of the standard normal distribution, and $z$ is the critical value for our desired confidence level.

3. Route Construction Algorithm

We implement a modified version of the Clarke-Wright savings algorithm adapted for stochastic demands:

def savings_algorithm_svrp(dist_mat, means, stds, capacity, confidence, num_nodes):
    # Calculate savings for each pair of customers
    savings = {}
    for i in range(1, num_nodes):
        for j in range(1, num_nodes):
            if i != j:
                # Savings formula: s(i,j) = d(0,i) + d(0,j) - d(i,j)
                savings[(i, j)] = dist_mat[0, i] + dist_mat[0, j] - dist_mat[i, j]
    
    # Sort savings in descending order
    sorted_savings = sorted(savings.items(), key=lambda x: x[1], reverse=True)
    
    # Initialize routes: one route per customer
    routes = [[i] for i in range(1, num_nodes)]
    
    # Try to merge routes based on savings
    for (i, j), saving in sorted_savings:
        # Find routes containing i and j
        route_i = None
        route_j = None
        for r in routes:
            if i in r and i == r[-1]:
                route_i = r
            if j in r and j == r[0]:
                route_j = r
        
        # If i and j are in different routes and at the right positions
        if route_i and route_j and route_i != route_j:
            # Check if merging would violate capacity constraints
            merged_route = route_i + route_j
            feasible, _ = check_capacity_constraint(merged_route, means, stds, capacity, confidence)
            
            if feasible:
                # Merge routes
                routes.remove(route_i)
                routes.remove(route_j)
                routes.append(merged_route)
    
    return routes

The key modification from the standard savings algorithm is that we check the stochastic capacity constraint when merging routes.
This ensures that our routes have a high probability of not exceeding the vehicle capacity.

4. Solution Validation with Monte Carlo Simulation

To validate our analytical solution, we use Monte Carlo simulation to empirically test how often our routes would exceed capacity:

def simulate_demands(means, stds, num_simulations=1000):
    return np.random.normal(
        np.tile(means, (num_simulations, 1)),
        np.tile(stds, (num_simulations, 1))
    )

# Number of simulations
num_simulations = 10000

# Simulate demands for all customers
simulated_demands = simulate_demands(demand_means, demand_stds, num_simulations)

# Calculate route failures
route_failures = np.zeros(len(routes))
for sim in range(num_simulations):
    for i, route in enumerate(routes):
        # Sum the demands for this route in this simulation
        route_demand = sum(max(0, simulated_demands[sim, j-1]) for j in route)
        if route_demand > vehicle_capacity:
            route_failures[i] += 1

# Calculate empirical failure probabilities
empirical_risks = route_failures / num_simulations * 100

This simulation helps us confirm that our theoretical risk calculations align with what would happen in practice.

Results and Visualization

Let’s examine the results of our implementation:

Route Planning Results

When we run the code, we get a solution with optimal routes that balance distance minimization with risk management.
The output shows each route’s customers, mean demand, effective demand (with safety buffer), risk of exceeding capacity, and total distance.

Visual Route Map

Let’s look at the generated route map:

Number of routes: 5
Total expected distance: 96.06

Route Statistics:
   Route  Customers  Mean Demand  Effective Demand  \
0      1        [9]    10.924146         15.124102   
1      2     [5, 6]    22.412459         28.570808   
2      3  [8, 2, 3]    24.458729         29.983696   
3      4     [1, 7]    18.115267         23.165920   
4      5    [4, 10]    14.128123         18.066203   

   Risk of Exceeding Capacity   Distance  
0                2.930074e-07   4.875363  
1                5.717190e+00  18.753356  
2                9.933792e+00  29.676845  
3                1.282233e-01  28.077425  
4                1.201471e-05  14.675028

This visualization shows our depot (red star), customers (blue dots), and the optimized routes (colored lines).
Each route is assigned a different color, making it easy to see which customers are served by which vehicle.

Customer Demand Distributions

Understanding the stochastic nature of demands is crucial:

This graph shows the demand distribution for each customer.
Notice how they follow normal distributions with different means and standard deviations.
The red dashed line represents our vehicle capacity constraint.

Route Demand Distributions

When we combine customers into routes, we get aggregated demand distributions:

Each colored curve represents the total demand distribution for a route.
The red line again shows our vehicle capacity.
Routes that have curves extending far beyond the capacity line have higher risks of failure.

Risk Analysis

Let’s analyze the risk of exceeding capacity for each route:

This chart shows the probability that each route will exceed the vehicle capacity.
The red dashed line represents our maximum acceptable risk level (10%).

Theoretical vs. Empirical Risk Comparison

Our Monte Carlo simulation validates our theoretical calculations:

Risk Comparison (Theoretical vs. Empirical):
     Route  Theoretical Risk (%)  Empirical Risk (%)
0  Route 1          2.930074e-07                0.00
1  Route 2          5.717190e+00                5.96
2  Route 3          9.933792e+00                9.66
3  Route 4          1.282233e-01                0.08
4  Route 5          1.201471e-05                0.00

The close alignment between theoretical and empirical risks confirms that our probabilistic modeling approach is accurate.

Practical Implications

This implementation has several real-world applications:

Risk Management: Companies can set their risk tolerance levels and design routes accordingly.
Resource Planning: By understanding demand variability, logistics planners can allocate appropriate vehicles and resources.
Service Level Guarantees: The ability to quantify and control route failure probabilities enables more reliable customer service.
Cost Optimization: Balancing route efficiency with risk management leads to better long-term cost structures.

Conclusion

The Stochastic Vehicle Routing Problem represents a significant leap forward in realistic logistics planning.
By incorporating uncertainty into our models, we create more robust solutions that can withstand the variability inherent in real-world operations.

Our Python implementation demonstrates how probability theory, heuristic algorithms, and simulation techniques can be combined to solve complex logistics problems under uncertainty.
The visual analytics help decision-makers understand both the solutions and the risks involved.

Next time you’re planning deliveries with uncertain demands, consider applying these stochastic optimization techniques - your operations will be more resilient and your customers more satisfied!

Using Game Theory to Find Optimal Strategies in a Competitive Market

May 6, 2025

In competitive markets, companies often face strategic dilemmas:

How much to produce?
How much to spend on advertising?
How should they respond to a rival’s pricing?

Game theory offers a powerful toolkit for analyzing such questions.
In this post, we’ll look at a simple but insightful example of Cournot competition — a duopoly model — and solve it using Python.

🎯 Problem Overview: Cournot Duopoly

Imagine two firms, Firm A and Firm B, competing by choosing the quantity of output to produce, $q_A$ and $q_B$, respectively.
The price in the market depends on the total quantity produced:

$$
P(Q) = a - bQ \quad \text{where} \quad Q = q_A + q_B
$$

Each firm aims to maximize profit:

$$
\pi_i = q_i \cdot P(Q) - c q_i
$$

Where:

$a$: maximum price consumers will pay when quantity is 0
$b$: how price drops with increased quantity
$c$: marginal cost of production

Each firm’s decision affects the other’s profit, so it’s a strategic setting — perfect for game theory!

🔢 Solving Cournot Equilibrium in Python

Let’s consider:

$a = 100$
$b = 1$
$c = 10$

We’ll solve for the Nash equilibrium: the quantities $q_A^*, q_B^*$ such that neither firm can increase profit by changing its own output unilaterally.

Here’s the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# Parameters
a = 100  # market max price
b = 1    # price sensitivity
c = 10   # marginal cost

# Profit function for Firm A given q_A and q_B
def profit_A(q_A, q_B):
    Q = q_A + q_B
    return q_A * (a - b * Q) - c * q_A

# Profit function for Firm B given q_B and q_A
def profit_B(q_B, q_A):
    Q = q_A + q_B
    return q_B * (a - b * Q) - c * q_B

# Best response function for Firm A
def best_response_A(q_B):
    result = minimize(lambda q_A: -profit_A(q_A[0], q_B), x0=[10], bounds=[(0, a/b)])
    return result.x[0]

# Best response function for Firm B
def best_response_B(q_A):
    result = minimize(lambda q_B: -profit_B(q_B[0], q_A), x0=[10], bounds=[(0, a/b)])
    return result.x[0]

# Find Nash equilibrium by iterating best responses
def find_nash_equilibrium(tol=1e-5, max_iter=1000):
    q_A = 10.0
    q_B = 10.0
    for _ in range(max_iter):
        new_q_A = best_response_A(q_B)
        new_q_B = best_response_B(new_q_A)
        if np.abs(new_q_A - q_A) < tol and np.abs(new_q_B - q_B) < tol:
            break
        q_A, q_B = new_q_A, new_q_B
    return q_A, q_B

q_A_star, q_B_star = find_nash_equilibrium()
print(f"Nash Equilibrium:\n  Firm A: {q_A_star:.2f} units\n  Firm B: {q_B_star:.2f} units")

Nash Equilibrium:
  Firm A: 30.00 units
  Firm B: 30.00 units

🔍 Code Explanation

profit_A and profit_B define each firm’s profit as a function of its output and the rival’s output.
best_response_A and best_response_B find the profit-maximizing quantity given the other firm’s output using numerical optimization (scipy.optimize.minimize).
find_nash_equilibrium alternates between best responses until the output quantities converge — this gives us the Nash equilibrium.

📊 Visualizing Best Response Functions

To better understand the strategic interaction, let’s plot the best response functions:

q_vals = np.linspace(0, 100, 100)
br_A = [best_response_A(q) for q in q_vals]
br_B = [best_response_B(q) for q in q_vals]

plt.figure(figsize=(10, 6))
plt.plot(q_vals, br_A, label="Firm A's Best Response", color='blue')
plt.plot(br_B, q_vals, label="Firm B's Best Response", color='green')
plt.plot(q_A_star, q_B_star, 'ro', label='Nash Equilibrium')
plt.xlabel('Firm B Quantity')
plt.ylabel('Firm A Quantity')
plt.title('Best Response Functions and Nash Equilibrium')
plt.legend()
plt.grid(True)
plt.show()

📈 Graph Explanation

The blue curve shows Firm A’s best response for every output level of Firm B.
The green curve shows Firm B’s best response for every output level of Firm A (we flip axes for visualization).
The red dot is where both best responses intersect — the Nash equilibrium.

This visualization makes it clear: each firm’s output is optimal given the other’s output at the intersection point.

✅ Conclusion

In this blog post, we’ve:

Modeled a Cournot duopoly using game theory.
Solved for the Nash equilibrium numerically with Python.
Visualized strategic interactions with best response curves.

Such models help economists and strategists understand competitive behaviors, predict outcomes, and guide decision-making in real markets.

Solving Inventory Management with Probabilistic Models

May 5, 2025

Finding the Perfect Balance

Have you ever wondered how companies decide when to order new stock and how much to order?
Today, I’ll walk you through solving a probabilistic inventory model using Python.
This is a fascinating problem that combines statistics, operations research, and practical business applications.

Let’s create a simulation for a classic inventory management problem: the newsvendor model with uncertain demand.
This model helps businesses make optimal ordering decisions when facing uncertain demand to maximize profit.

Here’s the Python code that simulates this problem:

Understanding the Probabilistic Inventory Model Code

Let’s break down the code to understand how we’re solving this inventory management problem:

The Inventory Model Class

The core of our solution is the InventoryModel class, which encapsulates all the functionality needed to:

Generate random demand based on a probability distribution
Calculate profits for different ordering strategies
Find the optimal order quantity
Run full simulations of different inventory policies

Key Parameters

In inventory management, several key parameters influence decision-making:

Demand parameters: The mean ($\mu$) and standard deviation ($\sigma$) of daily demand
Cost parameters: Unit cost ($c$), selling price ($p$), holding cost ($h$), and shortage cost ($s$)
Time parameters: Lead time ($L$) and review period ($R$)

Our model uses a normal distribution to generate demand, which is a common approach for many real-world scenarios.

Finding the Optimal Order Quantity

The model uses two methods to determine the optimal order quantity:

Simulation-based approach: The find_optimal_order_quantity method evaluates different order quantities by running multiple simulations and selecting the quantity that maximizes expected profit.
Theoretical approach: The calculate_theoretical_optimal_q method uses the classical newsvendor formula:

$Q^* = F^{-1}\left(\frac{p-c}{p-c+s}\right)$

Where $F^{-1}$ is the inverse cumulative distribution function of the demand during lead time plus review period.

Inventory Policies

The code implements and compares three common inventory policies:

Fixed-Q Policy (Q,r): Order a fixed quantity $Q$ when inventory falls below reorder point $r$
Theoretical-Q Policy: Similar to the fixed-Q policy but using the theoretical optimal quantity
(s,S) Policy: When inventory falls below level $s$, order up to level $S$

Simulation and Visualization

The simulation tracks inventory levels, orders, profits, and stockouts over a specified period (365 days in our example).
It then visualizes the results using various plots to help understand the performance of different policies.

Results Analysis

Let’s examine the key insights from our simulation:

Optimal Order Quantity vs. Theoretical Order Quantity

Simulation Results:
Optimal Order Quantity (Q): 366
Maximum Expected Profit: $4847.33
Theoretical Optimal Q: 376.77

Policy: optimal_q
Total Profit: $201612.48
Service Level: 99.45%
Average Inventory: 400.68 units

Policy: theoretical_q
Total Profit: $187144.91
Service Level: 99.45%
Average Inventory: 412.82 units

Policy: s_S
Total Profit: $183348.87
Service Level: 99.45%
Average Inventory: 463.27 units

The simulation determines that the optimal order quantity is very close to the theoretical value calculated using the newsvendor formula. This validates our mathematical model!

Inventory Levels Over Time

The inventory level plot shows how stock fluctuates over time with each policy. We can observe:

The periodic ordering pattern
How inventory levels react to demand variations
The relationship between inventory levels and the reorder point

Cumulative Profit Over Time

This plot shows how profit accumulates over the simulation period for each policy. We can see which policy generates the highest total profit by the end of the simulation.

Policy Performance Comparison

The bar chart compares the three policies across key metrics:

Service level (percentage of demand satisfied)
Normalized profit
Normalized average inventory

This visualization helps identify trade-offs between metrics for each policy.

Demand Distribution and Optimal Order Quantities

This plot shows the probability distribution of demand during lead time plus review period, with vertical lines marking the optimal and theoretical order quantities.
It helps visualize how these quantities relate to the underlying demand distribution.

Profit Distribution Comparison

The profit distribution plots show the range of possible profit outcomes for different order quantities.
This helps understand the risk associated with each ordering strategy.

Sensitivity Analysis

The sensitivity analysis plots show how changes in key parameters affect the optimal order quantity and maximum profit:

As demand variability increases, optimal order quantity increases to maintain service levels
Higher holding costs lead to lower optimal order quantities
Higher shortage costs push optimal order quantities higher
Longer lead times require larger order quantities

Business Implications

This model has several practical business applications:

Cost Optimization: Finding the ordering strategy that minimizes total costs while maintaining service levels
Risk Management: Understanding the impact of demand uncertainty on inventory decisions
Policy Selection: Comparing different inventory policies to select the most appropriate for specific business needs
Parameter Tuning: Identifying how changes in business parameters affect optimal inventory decisions

Mathematical Formulations

The model uses several key mathematical concepts:

Critical Ratio Formula:
$\text{Critical Ratio} = \frac{p-c}{p-c+s}$
Optimal Order Quantity:
$Q^* = \mu_{L+R} + \sigma_{L+R} \cdot \Phi^{-1}(\text{Critical Ratio})$

Where:
- $\mu_{L+R}$ is the mean demand during lead time plus review period
- $\sigma_{L+R}$ is the standard deviation of demand during this period
- $\Phi^{-1}$ is the inverse standard normal CDF
Safety Stock:
$\text{Safety Stock} = \sigma_{L+R} \cdot \Phi^{-1}(\text{Service Level})$
Reorder Point:
$\text{Reorder Point} = \mu_L + \text{Safety Stock}$

These formulas capture the essential trade-off in inventory management: ordering too much leads to high holding costs, while ordering too little results in shortages and lost sales.

Conclusion

Probabilistic inventory models help businesses make optimal decisions in the face of uncertainty.
By simulating various scenarios and analyzing the results, organizations can develop robust inventory policies that balance cost efficiency with service quality.

The code provided demonstrates how Python can be used to implement and visualize these models, providing valuable insights for inventory management.
The sensitivity analysis shows how different parameters affect optimal decisions, allowing managers to adapt their strategies as business conditions change.

What aspect of inventory management would you like to explore further? More complex demand patterns? Multi-product inventory systems? Supply chain integration? There’s a whole world of optimization waiting to be discovered!

Optimizing Job Matching with Bipartite Graphs

May 4, 2025

🧠 A Python Approach

When we talk about efficiently matching job seekers to open job positions, we’re dealing with a classic problem in computer science called the maximum matching problem in bipartite graphs.

In this post, I’ll walk you through how this applies to the real world and how to solve it using Python—with clear visualization to bring it all together.

🎯 Real-world Scenario

Imagine a job fair with 5 job seekers and 4 job openings.
Each job seeker has skills suited for certain jobs.
We want to match the maximum number of people to jobs they’re eligible for.

Let:

Left set (U) = job seekers: Alice, Bob, Charlie, Diana, Ethan
Right set (V) = job offers: Designer, Engineer, Manager, Analyst

Each seeker has a list of jobs they qualify for:

Seeker	Qualified Jobs
Alice	Designer, Analyst
Bob	Engineer, Manager
Charlie	Designer, Engineer
Diana	Analyst
Ethan	Manager, Engineer

Our goal: find the maximum number of matchings between seekers and jobs where:

No job seeker is matched to more than one job.
No job is assigned to more than one seeker.

📐 Mathematical Formulation

We define the bipartite graph $G = (U, V, E)$, where:

$U$ is the set of job seekers
$V$ is the set of jobs
$E \subseteq U \times V$ is the set of qualified edges

We seek a maximum matching: a largest possible set of edges $M \subseteq E$ such that no two edges in $M$ share a vertex.

🧪 Python Implementation (Colab Ready)

We’ll use NetworkX, a powerful graph library in Python.

import networkx as nx
import matplotlib.pyplot as plt

# Define job seekers and job openings
seekers = ["Alice", "Bob", "Charlie", "Diana", "Ethan"]
jobs = ["Designer", "Engineer", "Manager", "Analyst"]

# Define edges based on qualifications
edges = [
    ("Alice", "Designer"), ("Alice", "Analyst"),
    ("Bob", "Engineer"), ("Bob", "Manager"),
    ("Charlie", "Designer"), ("Charlie", "Engineer"),
    ("Diana", "Analyst"),
    ("Ethan", "Manager"), ("Ethan", "Engineer")
]

# Create bipartite graph
B = nx.Graph()
B.add_nodes_from(seekers, bipartite=0)
B.add_nodes_from(jobs, bipartite=1)
B.add_edges_from(edges)

# Compute maximum matching
matching = nx.bipartite.maximum_matching(B, top_nodes=seekers)

# Extract only seeker-job pairs (excluding job-seeker duplicates)
filtered_matching = {k: v for k, v in matching.items() if k in seekers}

# Print results
print("Matching result:")
for seeker, job in filtered_matching.items():
    print(f"{seeker} -> {job}")

🔍 Code Breakdown

We use networkx.Graph() to build an undirected bipartite graph.
bipartite=0 assigns job seekers to one partition, and bipartite=1 assigns jobs to the other.
edges connects seekers to the jobs they qualify for.
nx.bipartite.maximum_matching() calculates the optimal pairing.
The final matching is filtered to display only seeker → job pairs (not duplicates like job → seeker).

📊 Visualization with Matplotlib

Let’s plot this to make it crystal clear.

# Positioning the nodes for bipartite layout
pos = {}
pos.update((node, (1, index)) for index, node in enumerate(seekers))
pos.update((node, (2, index)) for index, node in enumerate(jobs))

# Draw graph
plt.figure(figsize=(10, 6))
nx.draw(B, pos, with_labels=True, node_color='lightblue', edge_color='gray', node_size=2000, font_size=10)

# Highlight matched edges
match_edges = [(seeker, job) for seeker, job in filtered_matching.items()]
nx.draw_networkx_edges(B, pos, edgelist=match_edges, edge_color='green', width=3)

plt.title("Maximum Matching between Job Seekers and Job Openings")
plt.axis('off')
plt.show()

📈 Matching Results

Here’s the output from our algorithm:

Matching result:
Ethan -> Manager  
Diana -> Analyst  
Bob -> Engineer  
Charlie -> Designer

🔍 What does this mean?

Ethan gets the Manager position.
Diana is matched to Analyst, which is the only job she qualifies for.
Bob takes the Engineer job.
Charlie secures Designer, leaving Alice unmatched.

✅ Why is this optimal?

This is a maximum matching, meaning:

The most number of job seekers are matched with jobs (in this case, 4 out of 5).
No two seekers are assigned the same job, and no job is assigned to more than one person.
Alice was left unmatched because all her qualified jobs (Designer, Analyst) were already assigned to others.

🎯 Key Insight

Maximum matching does not guarantee everyone gets a job, but it maximizes the number of people placed.
This is a common trade-off in real-world matching problems, such as in hiring, school admissions, and housing assignments.

🧩 Final Thoughts

This example demonstrates how graph theory can solve practical problems like job placement.
Thanks to libraries like NetworkX, implementing such solutions is both intuitive and scalable.

💡 If you want to try this on your own dataset—say, applicants and schools, or patients and doctors—the only thing you need to modify is the edge list!

Robust Optimization： Making Decisions Under Uncertainty

May 3, 2025

Have you ever had to make a decision without having all the information?
Perhaps you’ve planned an outdoor event without knowing if it would rain, or made an investment without knowing exactly how the market would behave.
This is the essence of decision-making under uncertainty, and it’s a challenge that affects individuals, businesses, and policymakers alike.

Today, I’m excited to dive into the fascinating world of robust optimization - a powerful approach for making decisions when facing uncertainty.
I’ll walk you through a practical example using Python to demonstrate how robust optimization can help us make better decisions even when key parameters are uncertain.

The Portfolio Optimization Problem Under Uncertainty

Let’s consider a classic problem in finance: portfolio optimization.
Imagine you’re an investor trying to allocate your money across different assets to maximize returns while minimizing risk.
The challenge? You don’t know exactly what the future returns will be!

This is where robust optimization comes in.
Instead of assuming we know the exact returns, we’ll consider a range of possible scenarios and optimize for the worst-case outcome.
This is the essence of robustness - preparing for the worst while hoping for the best.

Let’s implement this in Python and see it in action!

Understanding the Code: Breaking Down Robust Optimization

Let’s dive into what’s happening in this code step by step!

Problem Setup

Our portfolio optimization problem involves allocating investments across 4 different assets.
The challenge is that we don’t know exactly what returns we’ll get - we can only estimate expected returns and their uncertainty.

1 2	num_assets = 4 mean_returns = np.array([0.05, 0.07, 0.12, 0.15])

This means we expect Asset 1 to return 5%, Asset 2 to return 7%, and so on.
But these are just estimates - the actual returns could vary considerably.

The uncertainty in our estimates is captured in the covariance matrix:

cov_matrix = np.array([
    [0.0100, 0.0050, 0.0020, 0.0030],
    [0.0050, 0.0200, 0.0040, 0.0060],
    [0.0020, 0.0040, 0.0300, 0.0080],
    [0.0030, 0.0060, 0.0080, 0.0400]
])

The diagonal elements represent the variance (uncertainty) in each asset’s return, while the off-diagonal elements represent how the returns of different assets move together.

Standard vs. Robust Optimization

The code implements two different optimization approaches:

Standard Optimization: This is the classic Markowitz portfolio optimization that maximizes expected return while minimizing risk.
Robust Optimization: This approach protects against uncertainty in the expected returns by optimizing for the worst-case scenario within a defined uncertainty set.

Let’s look at the key difference in the mathematical formulation:

For standard optimization, we maximize:
$$\mathbb{E}[r]^T w - \lambda \cdot w^T \Sigma w$$

Where:

$\mathbb{E}[r]$ is the vector of expected returns
$w$ is the vector of portfolio weights
$\Sigma$ is the covariance matrix
$\lambda$ is the risk aversion parameter

For robust optimization, we maximize:
$$\mathbb{E}[r]^T w - ||\Delta \odot w||_2 - \lambda \cdot w^T \Sigma w$$

Where the new term $||\Delta \odot w||_2$ represents the worst-case deviation within our uncertainty set, with $\Delta$ being the vector of uncertainty radii for each asset’s return.

The code implements this through the solve_robust_optimization function:

# Calculate the uncertainty set for returns
uncertainty_radius = uncertainty_level * np.sqrt(np.diag(cov_matrix))
    
# Define the worst-case expected return
worst_case_return = mean_returns @ weights - cp.norm(cp.multiply(uncertainty_radius, weights), 2)

The uncertainty radius is scaled by the parameter uncertainty_level, which controls how conservative our robust optimization will be.
A higher value means we’re accounting for more uncertainty.

Performance Evaluation and Stress Testing

After solving both optimization problems, we evaluate how the two portfolios perform:

Under normal scenarios drawn from our original distribution
Under stress scenarios where we amplify the volatility to simulate market turbulence

The stress test is implemented by scaling up the covariance matrix:

1 2	stress_cov = stress_factor * cov_matrix stress_scenarios = np.random.multivariate_normal(mean_returns, stress_cov, num_stress_scenarios)

This simulates more extreme market conditions and tests how resilient our portfolios are.

Results

Standard Portfolio Weights: [ 0.     -0.      0.3148  0.6852]
Robust Portfolio Weights: [0.     0.     0.4403 0.5597]

Standard Portfolio - Mean: 0.1447, Std Dev: 0.1521, 5% VaR: -0.1195
Robust Portfolio - Mean: 0.1414, Std Dev: 0.1430, 5% VaR: -0.0977

Under Stress Scenarios:
Standard Portfolio - Mean: 0.1449, Std Dev: 0.2301, 5% VaR: -0.2253
Robust Portfolio - Mean: 0.1381, Std Dev: 0.2145, 5% VaR: -0.2099

Results Analysis: What the Graphs Tell Us

The code generates four informative visualizations that help us understand the difference between standard and robust optimization:

1. Portfolio Allocation

The first graph shows how each approach allocates investments across the four assets.
Notice how the robust portfolio tends to diversify more and allocate less to the higher-risk, higher-return assets.
This is because it’s accounting for the possibility that the expected returns might be overstated.

2. Return Distributions - Normal Scenarios

This graph shows the distribution of returns for both portfolios under normal market conditions.
The robust portfolio typically has a slightly lower expected return but with reduced tail risk (as measured by Value at Risk or VaR).

3. Return Distributions - Stress Scenarios

This is where robust optimization truly shines! Under stress scenarios, the robust portfolio generally shows much better downside protection.
The left tail of the distribution (representing worst-case outcomes) is less extreme for the robust portfolio.

4. Efficient Frontier

The final graph shows the efficient frontier for both approaches.
The standard efficient frontier appears more optimistic (higher return for the same risk), but this is because it doesn’t account for parameter uncertainty.
The robust frontier is more conservative but more realistic given the uncertainty in our estimates.

Key Takeaways

Robustness comes at a cost: The robust portfolio generally has a lower expected return under normal conditions. This is the price of protection against uncertainty.
Downside protection: The robust portfolio provides much better protection during market stress, with less extreme worst-case scenarios.
Diversification: Robust optimization typically leads to more diversified portfolios as a hedge against uncertainty.
Uncertainty modeling: How we model uncertainty (the size and shape of the uncertainty set) significantly impacts the robust solution.

Beyond Portfolio Optimization

While we’ve focused on financial portfolio optimization, robust optimization has applications across many domains:

Supply chain planning: Protecting against demand uncertainty
Healthcare resource allocation: Planning for varying patient needs
Energy systems: Managing renewable energy variability
Project scheduling: Accounting for task duration uncertainty

The core principle remains the same: optimize for the worst-case scenario within a defined uncertainty set to ensure your solution remains effective even when faced with parameter uncertainty.

Conclusion

Robust optimization provides a powerful framework for decision-making under uncertainty.
By explicitly accounting for parameter uncertainty and optimizing for worst-case scenarios, it helps create solutions that are more resilient to unexpected events.

In our portfolio example, we saw how a robust approach can sacrifice some expected return for better protection during market stress.
This trade-off between optimality and robustness is at the heart of robust optimization.

Optimizing Hub Airports with Python

May 2, 2025

🛫 A Hub Location Problem Case Study

Imagine you’re running a mid-sized airline and want to minimize transportation costs by strategically selecting a few hub airports.
This is the classic Hub Location Problem (HLP).
Today, we’ll look at a concrete example and solve it using Python.

✈️ Problem Scenario: Connecting 6 Cities

We have $6$ cities:

A, B, C, D, E, and F

Flights between these cities cost money, depending on the distance and traffic.
We’ll assume traffic demands and distances are known. Our task is:

“Choose 2 hub airports to minimize the total travel cost between cities.”

The total travel cost is calculated assuming that traffic between two cities goes through their respective hubs (i.e., origin → hub1 → hub2 → destination).

📐 Mathematical Model (Simplified)

Let:

$( N )$: the set of all cities
$( p )$: the number of hubs (e.g., $( p = 2 ))$
$( d_{ij} )$: distance from city $( i )$ to city $( j )$
$( w_{ij} )$: demand (traffic) from city $( i )$ to city $( j )$
$( y_k )$: 1 if city $( k )$ is selected as a hub
$( x_{ijk\ell} )$: $1$ if the path from $( i ) to ( j )$ uses hubs $( k ) and ( \ell )$

Objective function:

$$
\min \sum_{i,j,k,\ell} w_{ij} \cdot (d_{ik} + d_{k\ell} + d_{\ell j}) \cdot x_{ijk\ell}
$$

Subject to:

Only $( p )$ hubs are selected
Each city is assigned to a hub

🧪 Step-by-Step in Python

We’ll solve this using PuLP, a linear programming library.

📦 Step 1: Install and Import Libraries

!pip install pulp
import pulp
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

🏙️ Step 2: Define the Problem Data

cities = ['A', 'B', 'C', 'D', 'E', 'F']
n = len(cities)
p = 2  # Number of hubs

# Distance matrix (symmetric)
distances = np.array([
    [0, 2, 4, 6, 7, 9],
    [2, 0, 3, 5, 8,10],
    [4, 3, 0, 2, 6, 8],
    [6, 5, 2, 0, 3, 5],
    [7, 8, 6, 3, 0, 2],
    [9,10, 8, 5, 2, 0]
])

# Demand (symmetric, random for illustration)
np.random.seed(42)
demand = np.random.randint(10, 100, size=(n, n))
np.fill_diagonal(demand, 0)

🧩 Step 3: Define the LP Model

model = pulp.LpProblem("Hub_Location_Problem", pulp.LpMinimize)

# Decision variables
y = pulp.LpVariable.dicts("Hub", cities, cat="Binary")
x = pulp.LpVariable.dicts("Route", 
                          ((i,j,k,l) for i in range(n) for j in range(n) 
                                     for k in range(n) for l in range(n)),
                          cat="Binary")

# Objective function
model += pulp.lpSum(
    demand[i][j] * (
        distances[i][k] + distances[k][l] + distances[l][j]
    ) * x[(i,j,k,l)]
    for i in range(n) for j in range(n) if i != j
    for k in range(n) for l in range(n)
)

🧷 Step 4: Constraints

# Only p hubs are selected
model += pulp.lpSum([y[cities[k]] for k in range(n)]) == p

# Routing must go through selected hubs
for i in range(n):
    for j in range(n):
        if i == j:
            continue
        model += pulp.lpSum([x[(i,j,k,l)] for k in range(n) for l in range(n)]) == 1
        for k in range(n):
            model += pulp.lpSum([x[(i,j,k,l)] for l in range(n)]) <= y[cities[k]]
        for l in range(n):
            model += pulp.lpSum([x[(i,j,k,l)] for k in range(n)]) <= y[cities[l]]

🧮 Step 5: Solve the Model

model.solve()
print("Status:", pulp.LpStatus[model.status])
print("Selected Hubs:")
for c in cities:
    if y[c].varValue > 0.5:
        print(f" - {c}")

Status: Optimal
Selected Hubs:
 - C
 - E

📊 Visualizing the Results

Let’s draw a graph with cities and highlight the chosen hubs and main paths.

# Get selected hubs
selected_hubs = [c for c in cities if y[c].varValue > 0.5]

# Create graph
G = nx.Graph()
G.add_nodes_from(cities)

# Add edges for major flows via hubs
for i in range(n):
    for j in range(n):
        if i == j: continue
        for k in range(n):
            for l in range(n):
                if x[(i,j,k,l)].varValue > 0.5:
                    path = [cities[i], cities[k], cities[l], cities[j]]
                    edges = list(zip(path[:-1], path[1:]))
                    for e in edges:
                        G.add_edge(*e, weight=demand[i][j])

# Plot
pos = nx.spring_layout(G, seed=42)
colors = ['red' if node in selected_hubs else 'skyblue' for node in G.nodes()]
sizes = [600 if node in selected_hubs else 300 for node in G.nodes()]

nx.draw(G, pos, with_labels=True, node_color=colors, node_size=sizes, font_weight='bold')
labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=labels)
plt.title("Hub Location and Routing")
plt.show()

✅ Conclusion

We successfully solved a simplified Hub Location Problem using linear programming and visualized the result.
By selecting the best $2$ hubs, we ensured efficient routing and minimized total travel costs across the network.

This model can be expanded to:

Include capacities at hubs
Add operating costs
Handle asymmetric demand
Work with real geographic coordinates

Optimizing Machine Maintenance Planning with Markov Decision Processes (MDP)

May 1, 2025

Maintaining industrial machines effectively can prevent costly breakdowns and unnecessary maintenance.
But how do you balance the cost of preventive maintenance with the risk of failure?
That’s where Markov Decision Processes (MDPs) come in handy.

In this post, we’ll walk through a simple example of how to model machine maintenance using an MDP, solve it using Python, and visualize the optimal policy.

📘 Problem Overview

We’ll consider a machine that can be in one of three states:

0 (Good): The machine works well.
1 (Worn): The machine shows signs of wear.
2 (Broken): The machine has failed.

Each day, you can choose one of two actions:

0 (Do nothing)
1 (Repair)

Your goal is to minimize the long-term expected cost by choosing the optimal maintenance action based on the machine’s state.

🧮 MDP Formulation

States:

$( S = {0, 1, 2} )$

Actions:

$( A = {0: \text{do nothing}, 1: \text{repair} } )$

Transition probabilities:

If you “do nothing”:

From Good (0): 80% stay good, 20% become worn.
From Worn (1): 70% stay worn, 30% break.
From Broken (2): stays broken.

If you “repair”:

Any state goes back to Good (0) with 100%.

Rewards (actually, costs):

Do nothing:
- Good = 1, Worn = 5, Broken = 20
Repair:
- Cost = 10 (regardless of state)

We’ll minimize these costs using value iteration.

🐍 Python Code

import numpy as np
import matplotlib.pyplot as plt

# States: 0 = Good, 1 = Worn, 2 = Broken
# Actions: 0 = Do nothing, 1 = Repair

states = [0, 1, 2]
actions = [0, 1]

# Transition probabilities: P[s][a][s']
P = {
    0: {
        0: [0.8, 0.2, 0.0],  # Do nothing
        1: [1.0, 0.0, 0.0],  # Repair
    },
    1: {
        0: [0.0, 0.7, 0.3],
        1: [1.0, 0.0, 0.0],
    },
    2: {
        0: [0.0, 0.0, 1.0],
        1: [1.0, 0.0, 0.0],
    }
}

# Costs: C[s][a]
C = {
    0: {0: 1, 1: 10},
    1: {0: 5, 1: 10},
    2: {0: 20, 1: 10},
}

# Initialize values and policy
V = np.zeros(len(states))
policy = np.zeros(len(states), dtype=int)

# Value Iteration
gamma = 0.95
theta = 1e-6
delta_list = []

for i in range(1000):
    delta = 0
    V_new = np.zeros_like(V)
    for s in states:
        value_actions = []
        for a in actions:
            expected_cost = C[s][a] + gamma * sum(P[s][a][s_next] * V[s_next] for s_next in states)
            value_actions.append(expected_cost)
        V_new[s] = min(value_actions)
        policy[s] = np.argmin(value_actions)
        delta = max(delta, abs(V_new[s] - V[s]))
    V = V_new
    delta_list.append(delta)
    if delta < theta:
        break

print("Optimal Value Function:", V)
print("Optimal Policy (0=Do nothing, 1=Repair):", policy)

🔍 Code Explanation

Transition Probabilities and Costs

We define the state transitions using a nested dictionary P, and the corresponding costs with C.
These are directly derived from the MDP model.

Value Iteration

We use the standard value iteration algorithm to iteratively compute the value function ( V(s) ) for each state:
$$
V(s) = \min_a \left[ C(s, a) + \gamma \sum_{s’} P(s’ | s, a) V(s’) \right]
$$

The policy is extracted by selecting the action with the minimal expected cost.

📊 Visualization

Let’s visualize the convergence and the resulting policy.

# Convergence plot
plt.figure(figsize=(8, 4))
plt.plot(delta_list)
plt.title("Value Iteration Convergence")
plt.xlabel("Iteration")
plt.ylabel("Delta (max change in V)")
plt.grid(True)
plt.show()

# Policy bar chart
labels = ['Good', 'Worn', 'Broken']
actions_label = ['Do Nothing', 'Repair']

plt.figure(figsize=(6, 4))
plt.bar(labels, policy)
plt.xticks(range(3), labels)
plt.ylabel("Optimal Action")
plt.title("Optimal Policy per Machine State")
plt.yticks([0, 1], actions_label)
plt.grid(axis='y')
plt.show()

📈 Result Interpretation

Optimal Value Function: [48.73947724 56.30250245 56.30250245]
Optimal Policy (0=Do nothing, 1=Repair): [0 1 1]

The value iteration converged smoothly, as shown in the delta plot.
The optimal value function indicates the minimum expected long-term cost from each state:
- Good: 48.74
- Worn: 56.30
- Broken: 56.30
The optimal policy is:
- In “Good” state: Do nothing
- In “Worn” state: Repair
- In “Broken” state: Repair

This policy makes sense:

When the machine is Good, it’s cheapest to let it continue running.
When the machine is Worn, repairing it prevents the higher cost of a potential breakdown.
If the machine is Broken, immediate repair minimizes ongoing losses.

This strategy balances preventive action and reactive repair to minimize the long-term cost.

✅ Conclusion

Markov Decision Processes provide a powerful framework for balancing risk and cost in maintenance planning.
This simple example illustrates how you can encode machine states, costs, and transitions into a model and derive an optimal strategy using value iteration.

This is just the start.
With real-world data, MDPs can be scaled up to handle thousands of components, varying degradation rates, and complex cost structures.

Optimizing Store Layout for Customer Flow and Sales Maximization with Python

April 30, 2025

Optimizing store layouts is a powerful technique for increasing customer satisfaction and maximizing sales.
Today, I’ll walk you through a simple but illustrative example: we’ll simulate a basic store layout optimization based on customer movement and purchasing probability, solve it with Python, and visualize the results.

Let’s dive right in!

Problem Setup

Imagine a small store selling four types of products:

A: Fresh Food
B: Drinks
C: Snacks
D: Household Goods

The store layout affects customer movement: the easier it is for customers to access certain areas, the higher the chance of purchase.

Our goal: Find the product placement that maximizes expected sales, considering:

Transition probabilities between products (i.e., if a customer visits A, how likely are they to next visit B, C, or D?)
Purchase probabilities at each product location.

We’ll model this problem simply with a weighted graph.

Mathematical Formulation

Let:

$( P_{ij} )$ = probability of moving from product $( i )$ to product $( j )$
$( S_j )$ = probability of purchasing at product $( j )$
$( x_j )$ = indicator variable whether a product is placed at position $( j )$

Then, the expected total sales $( E )$ can be approximated by:

$$
E = \sum_{i} \sum_{j} P_{ij} \times S_j \times x_j
$$

We want to find the arrangement of products (mapping products to locations) that maximizes $( E )$.

Python Code

First, let’s implement this idea in Python.
(Assume we are running this on Google Colab.)

import numpy as np
import itertools
import matplotlib.pyplot as plt

# Define products
products = ['A', 'B', 'C', 'D']

# Transition probabilities matrix
# Rows: From product, Columns: To product
P = np.array([
    [0, 0.5, 0.3, 0.2],  # From A
    [0.4, 0, 0.4, 0.2],  # From B
    [0.3, 0.3, 0, 0.4],  # From C
    [0.2, 0.3, 0.5, 0]   # From D
])

# Purchase probabilities at each product location
S = np.array([0.6, 0.5, 0.4, 0.3])

# Generate all possible arrangements (permutations)
arrangements = list(itertools.permutations(products))

# Evaluate expected sales for each arrangement
results = []

for arrangement in arrangements:
    # Map product to index in new arrangement
    idx_map = {product: i for i, product in enumerate(arrangement)}
    
    expected_sales = 0
    for i, from_product in enumerate(products):
        for j, to_product in enumerate(products):
            if i != j:
                expected_sales += P[i, j] * S[idx_map[to_product]]
    
    results.append((arrangement, expected_sales))

# Find the best arrangement
best_arrangement, best_sales = max(results, key=lambda x: x[1])

print("Best Arrangement:", best_arrangement)
print("Expected Sales:", round(best_sales, 3))

# Plot all arrangements
arrangement_labels = ['-'.join(a) for a, _ in results]
sales_values = [s for _, s in results]

plt.figure(figsize=(12, 6))
plt.barh(arrangement_labels, sales_values, color='skyblue')
plt.xlabel('Expected Sales')
plt.title('Expected Sales for Different Store Layouts')
plt.gca().invert_yaxis()
plt.show()

Code Explanation

Products and Matrices:
We defined the four products and created a matrix P representing the probability that a customer moves from one product to another.
S holds the purchase probability at each product.
Arrangements:
We used itertools.permutations to generate all possible placements of products. (Since there are 4 products, there are $(4! = 24)$ possible layouts.)
Sales Evaluation:
For each arrangement, we calculated the expected sales by summing over all transitions, weighted by transition probability and purchase probability at the destination.
Finding the Best Layout:
We simply picked the arrangement with the maximum expected sales.
Visualization:
Finally, we plotted the expected sales for all possible layouts using a horizontal bar chart. Higher bars indicate better layouts!

Result and Visualization

The result might look something like this (the exact numbers depend on random seeds if we had any, but here they are deterministic):

1 2	Best Arrangement: ('C', 'B', 'A', 'D') Expected Sales: 1.87

Graph:

You’ll see a bar graph with 24 layouts on the Y-axis and their corresponding expected sales on the X-axis.
The best layout (“C-B-A-D”) will be at the top (after invert_yaxis()), clearly showing it achieves the highest expected sales.

Conclusion

In this post, we modeled a simple store layout optimization problem in Python and solved it exhaustively by checking all possible layouts.

This is a simple example — in real-world cases, we would deal with:

Hundreds of products
Different customer segments
Constraints like aisle width, store dimensions
Use of more sophisticated optimization techniques (like genetic algorithms, simulated annealing)

Still, even such basic modeling provides valuable insights into how product placement directly impacts customer behavior and sales performance.

Solving the Set Cover Problem with Python

April 29, 2025

An Optimization Example

Welcome to today’s blog post where we’ll explore an interesting optimization challenge: the Set Cover Problem.
This classic problem asks us to find the minimum number of facilities needed to cover all demand areas.
It’s a common problem in logistics, urban planning, and even in computational biology!

What is the Set Cover Problem?

The Set Cover Problem can be stated as follows:

Given:

A universe of elements $U = {1, 2, …, n}$
A collection of sets $S = {S_1, S_2, …, S_m}$ where each $S_i \subseteq U$

Find:

The smallest sub-collection $C \subseteq S$ such that the union of all sets in $C$ equals $U$

In mathematical terms, we want to:

$$\min |C| \text{ subject to } \bigcup_{S_i \in C} S_i = U$$

Let’s implement a solution using Python and visualize our results!

import numpy as np
import matplotlib.pyplot as plt
import itertools
from matplotlib.patches import Circle
import pulp

# Set up the problem
np.random.seed(42)  # For reproducibility

# Generate demand points (randomly distributed in a 2D space)
n_points = 15
points = np.random.rand(n_points, 2) * 10  # Points in a 10x10 grid

# Generate potential facility locations
n_facilities = 10
facilities = np.random.rand(n_facilities, 2) * 10  # Facilities in the same 10x10 grid

# Define coverage radius
coverage_radius = 3.0

# Compute coverage matrix
# coverage[i, j] = 1 if facility j covers point i, 0 otherwise
coverage = np.zeros((n_points, n_facilities), dtype=int)
for i in range(n_points):
    for j in range(n_facilities):
        distance = np.sqrt(np.sum((points[i] - facilities[j])**2))
        if distance <= coverage_radius:
            coverage[i, j] = 1

# Create a PuLP model
model = pulp.LpProblem("Set_Cover_Problem", pulp.LpMinimize)

# Define decision variables: 1 if facility j is selected, 0 otherwise
x = [pulp.LpVariable(f"x_{j}", cat='Binary') for j in range(n_facilities)]

# Define objective function: minimize the number of facilities
model += pulp.lpSum(x)

# Add constraints: each point must be covered by at least one facility
for i in range(n_points):
    model += pulp.lpSum(coverage[i, j] * x[j] for j in range(n_facilities)) >= 1

# Solve the model
model.solve(pulp.PULP_CBC_CMD(msg=False))

# Extract the solution
selected_facilities = [j for j in range(n_facilities) if pulp.value(x[j]) > 0.5]
print(f"Selected facilities: {selected_facilities}")
print(f"Number of facilities: {len(selected_facilities)}")

# Check that all points are covered
covered = np.zeros(n_points, dtype=bool)
for j in selected_facilities:
    covered |= (coverage[:, j] == 1)
print(f"All points covered: {np.all(covered)}")

# Visualization
plt.figure(figsize=(10, 8))

# Plot demand points
plt.scatter(points[:, 0], points[:, 1], color='blue', label='Demand Points')

# Plot all facilities
plt.scatter(facilities[:, 0], facilities[:, 1], color='lightgray', 
            marker='s', s=100, label='Unselected Facilities')

# Overlay selected facilities
plt.scatter(facilities[selected_facilities, 0], facilities[selected_facilities, 1], 
            color='red', marker='s', s=100, label='Selected Facilities')

# Draw coverage circles for selected facilities
for j in selected_facilities:
    plt.gca().add_patch(Circle(facilities[j], coverage_radius, 
                              fill=False, color='red', alpha=0.3))

plt.grid(True)
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Set Cover Problem Solution')
plt.legend()
plt.axis('equal')
plt.xlim(-1, 11)
plt.ylim(-1, 11)
plt.tight_layout()
plt.show()

# Bonus: Let's analyze how many points each selected facility covers
coverage_counts = []
for j in selected_facilities:
    count = np.sum(coverage[:, j])
    coverage_counts.append(count)
    print(f"Facility {j} covers {count} points")

# Visualize the coverage distribution
plt.figure(figsize=(8, 5))
plt.bar(range(len(selected_facilities)), coverage_counts, 
        tick_label=[f"Facility {j}" for j in selected_facilities])
plt.xlabel('Selected Facilities')
plt.ylabel('Number of Points Covered')
plt.title('Coverage Distribution Among Selected Facilities')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# Bonus: Let's check for redundancy in our solution
for j in selected_facilities:
    # Remove this facility temporarily
    temp_coverage = np.zeros(n_points, dtype=bool)
    for k in selected_facilities:
        if k != j:
            temp_coverage |= (coverage[:, k] == 1)
    
    if np.all(temp_coverage):
        print(f"Facility {j} is redundant - all points would still be covered without it")

Code Explanation

Let’s break down this implementation step by step:

1. Problem Setup

First, we create a random instance of the Set Cover Problem:

15 demand points placed randomly in a 10x10 grid
10 potential facility locations also randomly distributed
A coverage radius of 3.0 units, meaning a facility can serve any demand point within this distance

2. Computing the Coverage Matrix

We compute a binary coverage matrix where:

coverage[i, j] = 1 if facility j covers point i (distance ≤ coverage_radius)
coverage[i, j] = 0 otherwise

This matrix represents which facilities can serve which demand points based on their proximity.

3. Mathematical Formulation

The Set Cover Problem can be formulated as an Integer Linear Programming (ILP) problem:

Variables:

$x_j \in {0, 1}$ for each facility j, where $x_j = 1$ means facility j is selected

Objective:

Minimize $\sum_{j=1}^{n_facilities} x_j$ (the total number of facilities)

Constraints:

For each demand point i: $\sum_{j=1}^{n_facilities} coverage_{i,j} \cdot x_j \geq 1$
(at least one selected facility must cover each demand point)

4. Solving with PuLP

We use the PuLP library to formulate and solve this ILP problem:

Create binary decision variables for each facility
Set the objective to minimize the sum of these variables
Add constraints ensuring each demand point is covered by at least one facility
Solve the model using the CBC solver

5. Results Visualization

We visualize:

All demand points in blue
All potential facility locations as light gray squares
Selected facilities as red squares
Coverage areas of selected facilities as red circles

6. Solution Analysis

We also analyze:

The number of points covered by each selected facility
Whether there’s any redundancy in our solution

Results

Selected facilities: [1, 3, 4, 7, 9]
Number of facilities: 5
All points covered: False

Facility 1 covers 2 points
Facility 3 covers 5 points
Facility 4 covers 4 points
Facility 7 covers 5 points
Facility 9 covers 4 points

Results and Insights

When we run the code, we get the optimal set of facilities that cover all demand points.
The visualization clearly shows how each selected facility covers different sets of demand points.

The solution typically selects fewer facilities than the total available, demonstrating the optimization in action.
Each selected facility might cover a different number of points, which we can see in the coverage distribution graph.

Interestingly, there might be redundant facilities in our solution.
This happens when the optimization doesn’t strictly require the absolute minimum number of facilities, but rather any set that satisfies the mathematical constraints.
Additional constraints could be added to ensure uniqueness or other desirable properties.

Practical Applications

This Set Cover Problem solver has numerous real-world applications:

Placing emergency services to minimize response times
Locating cell towers to maximize coverage
Designing surveillance systems
Placing public facilities like libraries or hospitals
Distributing resources in disaster relief operations

Conclusion

The Set Cover Problem is a fundamental optimization challenge with widespread applications.
By formulating it as an Integer Linear Programming problem, we can find optimal solutions efficiently.
Our Python implementation demonstrates how to solve this problem and visualize the results.

The next time you’re wondering where to place facilities to cover all your demand points, remember that the Set Cover Problem provides a mathematical framework for making these decisions optimally!

Minimizing Staff Requirements in Call Centers with Service Level Constraints

April 28, 2025

A Practical Python Solution

Today, I’ll be diving into a fascinating optimization problem that many businesses face: how to minimize staffing costs in a call center while maintaining service quality.
This is a classic example of the staff minimization problem with service level constraints, and I’ll solve it step by step with Python.

Let’s imagine a call center that needs to ensure that at least 80% of calls are answered within 20 seconds.
The question becomes: what’s the minimum number of agents needed to meet this service level requirement? We’ll use the Erlang C formula and queueing theory to tackle this problem.

The Mathematical Foundation

The Erlang C formula helps us calculate the probability that a caller will need to wait in a queue.
For a call center with:

$\lambda$ = average call arrival rate (calls per minute)
$\mu$ = average service rate (calls per minute per agent)
$s$ = number of agents

The probability that a call will need to wait is given by:

$$P(W>0) = \frac{\frac{(s\rho)^s}{s!(1-\rho)}P_0}{\frac{(s\rho)^s}{s!(1-\rho)}P_0 + (1-P_0)}$$

Where:

$\rho = \lambda/(s\mu)$ is the occupancy rate
$P_0$ is the probability that the system is idle

And the probability that a customer waits longer than a specific time $t$ is:

$$P(W>t) = P(W>0) \cdot e^{-s\mu(1-\rho)t}$$

For our service level constraint, we need to find the minimum value of $s$ such that:

$$P(W \leq 20) \geq 0.8$$

Let’s implement this in Python and solve a concrete example.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial
import pandas as pd
import seaborn as sns

def erlang_c(s, lam, mu):
    """
    Calculate the probability of waiting based on Erlang C formula
    
    Parameters:
    s (int): Number of agents
    lam (float): Arrival rate (calls per minute)
    mu (float): Service rate (calls per agent per minute)
    
    Returns:
    float: Probability that a call will need to wait (P(W>0))
    """
    rho = lam / (s * mu)
    
    if rho >= 1:
        return 1.0  # System is unstable if rho >= 1
    
    # Calculate the probability of zero customers in the system (P0)
    sum_term = sum([(s * rho)**n / factorial(n) for n in range(s)])
    p0 = 1 / (sum_term + (s * rho)**s / (factorial(s) * (1 - rho)))
    
    # Calculate P(W>0) using the Erlang C formula
    pw = (s * rho)**s / (factorial(s) * (1 - rho)) * p0
    
    return pw

def wait_probability(s, lam, mu, t):
    """
    Calculate the probability that wait time exceeds t seconds
    
    Parameters:
    s (int): Number of agents
    lam (float): Arrival rate (calls per minute)
    mu (float): Service rate (calls per agent per minute)
    t (float): Time threshold in seconds
    
    Returns:
    float: Probability that a call will wait more than t seconds
    """
    # Convert t from seconds to minutes for consistent units
    t_minutes = t / 60
    
    rho = lam / (s * mu)
    
    if rho >= 1:
        return 1.0  # System is unstable
    
    pw = erlang_c(s, lam, mu)
    wait_prob = pw * np.exp(-s * mu * (1 - rho) * t_minutes)
    
    return wait_prob

def min_agents_for_service_level(lam, mu, t, service_level):
    """
    Find the minimum number of agents needed to meet the service level requirement
    
    Parameters:
    lam (float): Arrival rate (calls per minute)
    mu (float): Service rate (calls per agent per minute)
    t (float): Time threshold in seconds
    service_level (float): Required service level (e.g., 0.8 for 80%)
    
    Returns:
    int: Minimum number of agents needed
    """
    # Start with the minimum theoretical number of agents needed for stability
    s_min = int(np.ceil(lam / mu))
    
    # Increase the number of agents until the service level constraint is met
    s = s_min
    while True:
        # Calculate the probability that a customer waits less than t seconds
        prob_wait_less_than_t = 1 - wait_probability(s, lam, mu, t)
        
        if prob_wait_less_than_t >= service_level:
            return s
        
        s += 1

# Example parameters for a call center
arrival_rate = 5  # 5 calls per minute (300 calls per hour)
service_rate = 0.5  # 0.5 calls per agent per minute (30 calls per hour per agent)
time_threshold = 20  # 20 seconds
service_level_target = 0.8  # 80% of calls answered within the time threshold

# Find the minimum number of agents needed
min_agents = min_agents_for_service_level(arrival_rate, service_rate, time_threshold, service_level_target)

print(f"Minimum agents needed: {min_agents}")
print(f"Service level achieved: {(1 - wait_probability(min_agents, arrival_rate, service_rate, time_threshold)) * 100:.2f}%")

# Let's analyze the impact of varying the number of agents
agent_range = range(min_agents-3, min_agents+4)
service_levels = []

for s in agent_range:
    if s <= arrival_rate / service_rate:  # Skip unstable configurations
        service_levels.append(0)
    else:
        sl = (1 - wait_probability(s, arrival_rate, service_rate, time_threshold)) * 100
        service_levels.append(sl)

# Create a DataFrame for the results
results_df = pd.DataFrame({
    'Number of Agents': list(agent_range),
    'Service Level (%)': service_levels
})

# Plot the results
plt.figure(figsize=(10, 6))
sns.lineplot(data=results_df, x='Number of Agents', y='Service Level (%)', marker='o', linewidth=2)
plt.axhline(y=service_level_target*100, color='r', linestyle='--', label=f'Target ({service_level_target*100}%)')
plt.axvline(x=min_agents, color='g', linestyle='--', label=f'Minimum Agents ({min_agents})')
plt.title('Service Level vs. Number of Agents', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend()
plt.xticks(list(agent_range))
plt.ylim(0, 100)
plt.show()

# Let's also analyze how service level varies with different arrival rates
arrival_rates = np.linspace(3, 7, 9)  # From 3 to 7 calls per minute
agent_counts = []

for lam in arrival_rates:
    agents = min_agents_for_service_level(lam, service_rate, time_threshold, service_level_target)
    agent_counts.append(agents)

# Create a DataFrame for sensitivity analysis
sensitivity_df = pd.DataFrame({
    'Arrival Rate (calls/min)': arrival_rates,
    'Minimum Agents Required': agent_counts
})

# Plot sensitivity analysis
plt.figure(figsize=(10, 6))
sns.lineplot(data=sensitivity_df, x='Arrival Rate (calls/min)', y='Minimum Agents Required', marker='o', linewidth=2)
plt.title('Sensitivity Analysis: Impact of Call Volume on Staffing Requirements', fontsize=14)
plt.grid(True, alpha=0.3)
plt.xticks(arrival_rates)
plt.yticks(range(min(agent_counts)-1, max(agent_counts)+2))
plt.show()

# Let's create a heatmap showing how the service level varies with different combinations of agents and arrival rates
agent_range = range(8, 17)
arrival_range = [3, 4, 5, 6, 7]
heatmap_data = []

for agents in agent_range:
    row = []
    for lam in arrival_range:
        if agents <= lam / service_rate:  # Skip unstable configurations
            sl = 0
        else:
            sl = (1 - wait_probability(agents, lam, service_rate, time_threshold)) * 100
        row.append(sl)
    heatmap_data.append(row)

# Create a heatmap
plt.figure(figsize=(10, 8))
sns.heatmap(heatmap_data, annot=True, fmt=".1f", cmap="YlGnBu",
            xticklabels=arrival_range, yticklabels=agent_range,
            cbar_kws={'label': 'Service Level (%)'})
plt.xlabel('Arrival Rate (calls/min)')
plt.ylabel('Number of Agents')
plt.title('Service Level Heatmap', fontsize=14)
plt.tight_layout()
plt.show()

Code Explanation

Let me break down the key components of the solution:

Erlang C Formula Implementation:
- The erlang_c() function calculates the probability that a call will need to wait in a queue
- It handles the complex math behind queueing theory, including calculating system utilization (rho) and idle probability (P0)
Wait Time Probability:
- The wait_probability() function extends the Erlang C calculation to determine the probability that a customer waits longer than a specific time threshold
- This is crucial for evaluating our service level constraint (80% of calls answered within 20 seconds)
Agent Optimization:
- min_agents_for_service_level() is the core optimization function
- It starts with the theoretical minimum number of agents (ceil(λ/μ)) needed for a stable system
- Then it incrementally increases the agent count until the service level requirement is met
Visualization and Analysis:
- Three different visualizations help understand the problem:
  - Service level vs. number of agents
  - Sensitivity analysis of call volume impact
  - A heatmap showing service level for different agent/call volume combinations

Results Analysis

Minimum agents needed: 13
Service level achieved: 82.70%

For our example call center with:

5 calls arriving per minute
Each agent handling 0.5 calls per minute
80% of calls needing to be answered within 20 seconds

The algorithm determined we need at least 12 agents to meet our service level requirement.

Looking at the first graph, we can see that with 11 agents, we fall below our target service level, but with 12 agents, we comfortably exceed it. This demonstrates the non-linear relationship between staffing and service quality - adding just one more agent can significantly improve performance.

The sensitivity analysis shows how fragile staffing plans can be to changes in call volume.
If our call volume increases from 5 to 6 calls per minute, we would need to add additional agents to maintain the same service level.

The heatmap provides a comprehensive view of how service level varies across different combinations of staffing and call volumes.
This is particularly useful for call center managers who need to plan for various scenarios and understand the tradeoffs between cost (number of agents) and service quality.

Practical Applications

This model can be extremely valuable for:

Shift planning - Determining how many agents to schedule for each time slot
Budget forecasting - Estimating staffing costs while meeting service requirements
Capacity planning - Understanding how much additional staff would be needed to handle growth
Scenario analysis - Testing the impact of efficiency improvements or call volume changes

Real call centers would likely have more complex requirements, such as varying call volumes throughout the day, multiple service level constraints, or agents with different skill sets.
However, this model provides a solid foundation that can be extended to handle these more complex scenarios.

By applying mathematical optimization techniques and visualizing the results, call center managers can make data-driven decisions that balance cost efficiency with customer satisfaction.