Understanding Network Importance with PageRank

May 10, 2025

A Python Implementation

Have you ever wondered how Google decides which websites are more important than others?
Today, we’re diving into one of the most revolutionary algorithms in search engine history: PageRank.
Created by Larry Page and Sergey Brin, the founders of Google, this algorithm fundamentally changed how we navigate the internet.

What is PageRank?

PageRank is an algorithm that measures the importance of nodes in a network based on the structure of incoming links.
In simple terms, a page is important if other important pages link to it.
This creates a recursive definition that can be solved mathematically.

The core equation behind PageRank can be expressed as:

$$PR(A) = (1-d) + d \sum_{i=1}^n \frac{PR(T_i)}{C(T_i)}$$

Where:

$PR(A)$ is the PageRank of page A
$d$ is a damping factor (typically 0.85)
$T_i$ are pages linking to page A
$C(T_i)$ is the number of outbound links from page $T_i$

Let’s implement this algorithm in Python to analyze a simple network!

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.colors import LinearSegmentedColormap

# Set the style for our plots
plt.style.use('fivethirtyeight')
sns.set_palette("viridis")

def create_sample_network():
    """
    Create a sample directed network representing a citation network
    between academic papers or web pages.
    """
    # Create a directed graph
    G = nx.DiGraph()
    
    # Add nodes (representing web pages or academic papers)
    nodes = ["A", "B", "C", "D", "E", "F", "G", "H"]
    G.add_nodes_from(nodes)
    
    # Add edges (representing links between pages)
    edges = [
        ("A", "B"), ("A", "C"), ("A", "D"),
        ("B", "C"), ("B", "E"),
        ("C", "A"), ("C", "F"),
        ("D", "B"), ("D", "G"),
        ("E", "D"), ("E", "H"),
        ("F", "G"), ("F", "H"),
        ("G", "E"),
        ("H", "G")
    ]
    G.add_edges_from(edges)
    
    return G

def calculate_pagerank_matrix(G, alpha=0.85, max_iter=100, tol=1e-6):
    """
    Calculate PageRank using the power iteration method.
    
    Parameters:
    - G: NetworkX DiGraph
    - alpha: damping factor (default 0.85)
    - max_iter: maximum number of iterations (default 100)
    - tol: convergence tolerance (default 1e-6)
    
    Returns:
    - PageRank values dictionary
    - Number of iterations performed
    - Convergence history
    """
    n = len(G)
    
    # Create adjacency matrix
    A = nx.to_numpy_array(G, nodelist=list(G.nodes())).T
    
    # Create out-degree matrix
    out_degrees = np.array([max(G.out_degree(node), 1) for node in G.nodes()])
    D_inv = np.diag(1.0 / out_degrees)
    
    # Create transition probability matrix
    M = A @ D_inv
    
    # Initialize PageRank vector
    pr = np.ones(n) / n
    
    # For tracking convergence
    convergence_history = [np.copy(pr)]
    
    # Power iteration method
    for iteration in range(max_iter):
        prev_pr = np.copy(pr)
        
        # PageRank formula: (1-alpha)/n + alpha * M * pr
        pr = (1 - alpha) / n + alpha * M @ pr
        
        # Normalize to ensure sum equals 1
        pr = pr / np.sum(pr)
        
        # Store for convergence history
        convergence_history.append(np.copy(pr))
        
        # Check for convergence
        if np.linalg.norm(pr - prev_pr, 1) < tol:
            break
    
    # Create a dictionary with node names and their PageRank values
    pagerank = {node: rank for node, rank in zip(G.nodes(), pr)}
    
    return pagerank, iteration + 1, convergence_history

def visualize_network_with_pagerank(G, pagerank):
    """
    Visualize the network with node sizes proportional to PageRank values.
    """
    plt.figure(figsize=(12, 10))
    
    # Create a custom colormap from light to dark blue
    colors = ["#9ecae1", "#6baed6", "#4292c6", "#2171b5", "#084594"]
    custom_cmap = LinearSegmentedColormap.from_list("custom_blues", colors)
    
    # Node positions using spring layout
    pos = nx.spring_layout(G, seed=42)
    
    # Scale PageRank values for visualization
    node_sizes = [pagerank[node] * 10000 for node in G.nodes()]
    node_colors = [pagerank[node] for node in G.nodes()]
    
    # Draw nodes
    nodes = nx.draw_networkx_nodes(
        G, pos, 
        node_size=node_sizes,
        node_color=node_colors,
        cmap=custom_cmap,
        alpha=0.9
    )
    
    # Draw edges with arrows
    nx.draw_networkx_edges(
        G, pos, 
        width=1.5, 
        alpha=0.7, 
        edge_color='gray',
        arrowsize=15,
        connectionstyle='arc3,rad=0.1'
    )
    
    # Draw labels
    nx.draw_networkx_labels(
        G, pos, 
        font_size=12, 
        font_weight='bold',
        font_family='sans-serif'
    )
    
    # Add a colorbar
    plt.colorbar(nodes, label='PageRank Score', shrink=0.8)
    
    plt.title("Network Visualization with PageRank Importance", fontsize=18)
    plt.axis('off')
    plt.tight_layout()
    plt.savefig('network_pagerank.png', dpi=300, bbox_inches='tight')
    plt.show()
    
def visualize_pagerank_distribution(pagerank):
    """
    Create a bar chart showing the PageRank distribution.
    """
    plt.figure(figsize=(12, 6))
    
    # Sort PageRank values
    sorted_pagerank = {k: v for k, v in sorted(pagerank.items(), key=lambda item: item[1], reverse=True)}
    
    # Plot bar chart
    bars = plt.bar(sorted_pagerank.keys(), sorted_pagerank.values(), color=sns.color_palette("viridis", len(pagerank)))
    
    # Add value labels on top of bars
    for bar in bars:
        height = bar.get_height()
        plt.text(
            bar.get_x() + bar.get_width()/2., 
            height * 1.01,
            f'{height:.4f}',
            ha='center', 
            va='bottom', 
            fontsize=10,
            rotation=0
        )
    
    plt.title('PageRank Score Distribution', fontsize=18)
    plt.xlabel('Node', fontsize=14)
    plt.ylabel('PageRank Score', fontsize=14)
    plt.grid(axis='y', linestyle='--', alpha=0.7)
    plt.tight_layout()
    plt.savefig('pagerank_distribution.png', dpi=300, bbox_inches='tight')
    plt.show()

def plot_convergence(convergence_history, node_names):
    """
    Plot the convergence of PageRank values over iterations.
    """
    plt.figure(figsize=(12, 6))
    
    # Convert list of arrays to 2D array for easier indexing
    convergence_data = np.array(convergence_history)
    
    # Plot each node's PageRank convergence
    for i, node in enumerate(node_names):
        plt.plot(convergence_data[:, i], marker='o', markersize=4, label=f'Node {node}')
    
    plt.title('PageRank Convergence Over Iterations', fontsize=18)
    plt.xlabel('Iteration', fontsize=14)
    plt.ylabel('PageRank Value', fontsize=14)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.legend(loc='center right')
    plt.tight_layout()
    plt.savefig('pagerank_convergence.png', dpi=300, bbox_inches='tight')
    plt.show()

def comparative_analysis(G):
    """
    Compare our PageRank implementation with NetworkX's built-in PageRank.
    """
    # Calculate PageRank using our implementation
    our_pagerank, iterations, _ = calculate_pagerank_matrix(G)
    
    # Calculate PageRank using NetworkX's implementation
    nx_pagerank = nx.pagerank(G, alpha=0.85)
    
    # Create a DataFrame for comparison
    comparison_data = {
        'Node': list(our_pagerank.keys()),
        'Our Implementation': list(our_pagerank.values()),
        'NetworkX Implementation': [nx_pagerank[node] for node in our_pagerank.keys()]
    }
    
    # Plot comparison
    plt.figure(figsize=(12, 6))
    bar_width = 0.35
    x = np.arange(len(comparison_data['Node']))
    
    # Plot bars
    plt.bar(x - bar_width/2, comparison_data['Our Implementation'], bar_width, label='Our Implementation')
    plt.bar(x + bar_width/2, comparison_data['NetworkX Implementation'], bar_width, label='NetworkX Implementation')
    
    # Add details
    plt.title('PageRank Implementation Comparison', fontsize=18)
    plt.xlabel('Node', fontsize=14)
    plt.ylabel('PageRank Score', fontsize=14)
    plt.xticks(x, comparison_data['Node'])
    plt.legend()
    plt.grid(axis='y', linestyle='--', alpha=0.7)
    plt.tight_layout()
    plt.savefig('pagerank_comparison.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # Calculate absolute difference
    diff = np.abs(np.array(list(our_pagerank.values())) - np.array(list(nx_pagerank.values())))
    max_diff = np.max(diff)
    
    print(f"Maximum absolute difference between implementations: {max_diff:.10f}")
    print(f"Number of iterations until convergence: {iterations}")
    
    return our_pagerank, nx_pagerank

# Main execution
def main():
    # Create our sample network
    G = create_sample_network()
    
    # Calculate PageRank
    pagerank, iterations, convergence_history = calculate_pagerank_matrix(G)
    
    # Print results
    print("PageRank Results:")
    for node, rank in sorted(pagerank.items(), key=lambda x: x[1], reverse=True):
        print(f"Node {node}: {rank:.6f}")
    
    print(f"\nNumber of iterations until convergence: {iterations}")
    
    # Visualize the network with PageRank values
    visualize_network_with_pagerank(G, pagerank)
    
    # Visualize PageRank distribution
    visualize_pagerank_distribution(pagerank)
    
    # Plot convergence
    plot_convergence(convergence_history, list(G.nodes()))
    
    # Compare with NetworkX implementation
    our_pagerank, nx_pagerank = comparative_analysis(G)

if __name__ == "__main__":
    main()```


## Code Explanation

Let's break down the key components of the code:

### Network Creation

We start by creating a sample directed network with 8 nodes (A through H), where each node could represent a web page or an academic paper. 
The edges represent links or citations between these nodes.

```python
def create_sample_network():
    G = nx.DiGraph()
    
    # Add nodes (representing web pages or academic papers)
    nodes = ["A", "B", "C", "D", "E", "F", "G", "H"]
    G.add_nodes_from(nodes)
    
    # Add edges (representing links between pages)
    edges = [
        ("A", "B"), ("A", "C"), ("A", "D"),
        ("B", "C"), ("B", "E"),
        # ... other connections
    ]
    G.add_edges_from(edges)
    
    return G

This creates a directed graph where, for example, node A points to nodes B, C, and D.
This could represent a website A linking to websites B, C, and D.

PageRank Implementation

The core of our implementation is the calculate_pagerank_matrix function:

1	def calculate_pagerank_matrix(G, alpha=0.85, max_iter=100, tol=1e-6):

Here’s how it works:

We create an adjacency matrix A from the graph
We calculate the out-degree matrix D_inv, which contains the reciprocal of the number of outgoing links from each node
We compute the transition probability matrix M = A @ D_inv
We initialize the PageRank vector with equal probabilities
We iterate until convergence using the power iteration method

The main PageRank update equation in code:

1 2	# PageRank formula: (1-alpha)/n + alpha * M * pr pr = (1 - alpha) / n + alpha * M @ pr

This corresponds directly to the mathematical formula:

$$PR(A) = (1-d) + d \sum_{i=1}^n \frac{PR(T_i)}{C(T_i)}$$

The damping factor alpha (typically 0.85) represents the probability that a random surfer will follow a link rather than jump to a random page.

Visualization Functions

We’ve included several visualization functions to help understand the results:

visualize_network_with_pagerank: Displays the network with node sizes proportional to their PageRank scores
visualize_pagerank_distribution: Creates a bar chart showing the distribution of PageRank values
plot_convergence: Shows how PageRank values converge over iterations
comparative_analysis: Compares our implementation with NetworkX’s built-in PageRank function

Running the Algorithm

Let’s execute the code and analyze the results:

PageRank Results:
Node E: 0.243118
Node G: 0.218767
Node H: 0.142874
Node D: 0.135941
Node B: 0.090391
Node C: 0.071032
Node A: 0.048939
Node F: 0.048939

Number of iterations until convergence: 63

Maximum absolute difference between implementations: 0.0000020290
Number of iterations until convergence: 63

Analyzing the Results

Let’s analyze the output and visualize what our PageRank algorithm has discovered about this network.

Network Visualization

In this visualization, the size and color intensity of each node corresponds to its PageRank value.
Nodes E and G have the highest PageRank values, indicating they are the most important nodes in our network.
This makes sense because:

PageRank Distribution

The bar chart clearly shows the distribution of PageRank values across all nodes.
We can see that:

Node E has the highest PageRank value (0.2431)
Node G follows closely (0.2188)
Node F has the lowest PageRank (0.0489)

This distribution reflects the structure of our network and how information or importance flows through it.

Convergence Analysis

The convergence plot shows how the PageRank values stabilize over iterations.
Most nodes reach a stable value after about 15-20 iterations.

This demonstrates the iterative nature of the PageRank algorithm and how it gradually refines its estimates of node importance.

Implementation Comparison

Our implementation produces identical results to NetworkX’s built-in PageRank function, with a maximum absolute difference of practically zero.
This validates our implementation and confirms that it correctly applies the PageRank algorithm principles.

Why Does This Matter?

PageRank has applications far beyond web search:

Academic Citation Networks: Identifying influential papers
Social Network Analysis: Finding key influencers
Biological Networks: Understanding protein-protein interactions
Infrastructure Analysis: Identifying critical nodes in transportation networks
Recommendation Systems: Ranking items based on user preferences

Conclusion

The PageRank algorithm is a powerful tool for analyzing network structures and determining node importance.
Our implementation demonstrates how relatively simple matrix operations can reveal complex network dynamics.

In our example network, nodes E and G emerged as the most important, while A and F were less significant.
These insights could help prioritize resources, identify influential actors, or understand information flow within the network.

The next time you search for something on Google, remember that a variation of this algorithm is working behind the scenes to rank the pages you see!

Minimizing Cost While Ensuring Service Level

May 9, 2025

📦 Optimal Inventory Management with Python

Inventory management is a delicate balancing act.
Too much stock ties up capital and raises storage costs; too little risks lost sales and dissatisfied customers.
In this post, we’ll walk through a real-world example using Python to find the optimal inventory level that minimizes cost while maintaining a specific service level.

We’ll use a classical inventory model that considers:

Holding cost
Stockout cost (penalty)
Service level target (e.g. 95%)

Let’s dive in.

📘 Problem Setup

Imagine you’re managing a warehouse for a retail company.
A particular product has:

Daily demand: Normally distributed with a mean of 100 units and standard deviation of 20.
Lead time: 5 days
Holding cost: $2 per unit per day
Stockout cost: $50 per unit short
Target service level: 95%

We want to find the optimal reorder point $R$ and order quantity $Q$ such that we:

Minimize the total expected cost
Maintain a 95% service level (i.e., only 5% of demand during lead time results in stockouts)

🔢 Mathematical Model

Let’s define:

$\mu$: Average daily demand
$\sigma$: Standard deviation of daily demand
$L$: Lead time (days)
$h$: Holding cost per unit per day
$p$: Stockout penalty per unit
$SL$: Desired service level

Lead time demand is normally distributed:

Mean: $\mu_L = \mu \cdot L$
Std Dev: $\sigma_L = \sigma \cdot \sqrt{L}$

To achieve the desired service level $SL$, the reorder point $R$ is calculated as:

$$
R = \mu_L + z \cdot \sigma_L
$$

Where $z$ is the z-score corresponding to $SL$.

🧪 Python Implementation

Let’s translate all this into Python:

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# Parameters
mu = 100            # average daily demand
sigma = 20          # standard deviation of daily demand
L = 5               # lead time in days
h = 2               # holding cost per unit per day
p = 50              # stockout penalty per unit
service_level = 0.95

# Lead time demand distribution
mu_L = mu * L
sigma_L = sigma * np.sqrt(L)

# Find z-score for target service level
z = stats.norm.ppf(service_level)

# Reorder point R
R = mu_L + z * sigma_L

# Safety stock
safety_stock = z * sigma_L

# Expected shortage (standard normal loss function)
L_z = stats.norm.pdf(z) - z * (1 - service_level)

# Expected units short per cycle
E_short = sigma_L * L_z

# Total expected cost function over order quantity Q
def total_cost(Q):
    cycle_stock = Q / 2
    holding_cost = h * (cycle_stock + safety_stock)
    stockout_cost = p * E_short * (mu / Q)
    return holding_cost + stockout_cost

# Search for optimal order quantity Q
Q_values = np.arange(50, 1000, 10)
costs = [total_cost(Q) for Q in Q_values]
optimal_index = np.argmin(costs)
optimal_Q = Q_values[optimal_index]
min_cost = costs[optimal_index]

print(f"Reorder Point (R): {R:.2f}")
print(f"Safety Stock: {safety_stock:.2f}")
print(f"Optimal Order Quantity (Q): {optimal_Q}")
print(f"Minimum Total Cost: ${min_cost:.2f}")

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(Q_values, costs, label="Total Cost", color='blue')
plt.axvline(optimal_Q, color='red', linestyle='--', label=f"Optimal Q = {optimal_Q}")
plt.title("Total Expected Cost vs. Order Quantity")
plt.xlabel("Order Quantity (Q)")
plt.ylabel("Total Expected Cost ($)")
plt.legend()
plt.grid(True)
plt.show()

🧠 Code Breakdown

Input parameters define the demand distribution, lead time, and costs.
We calculate the lead time demand and z-score for the desired service level.
The reorder point $R$ ensures we cover average demand plus safety stock.
L_z computes the expected shortage per cycle using the standard normal loss function.
The total_cost() function combines:
- Cycle stock holding cost: $\frac{Q}{2} \cdot h$
- Safety stock holding cost: $SS \cdot h$
- Expected stockout cost per cycle
Finally, we search for the order quantity $Q$ that minimizes total cost using a simple grid search.

📊 Results Visualization

Reorder Point (R): 573.56
Safety Stock: 73.56
Optimal Order Quantity (Q): 70
Minimum Total Cost: $283.86

The graph above shows the relationship between order quantity $Q$ and total cost.
The red dashed line indicates the optimal order quantity, where the cost is minimized.

This visualization helps decision-makers see the cost sensitivity around the optimal point.

✅ Conclusion

With just a few lines of Python, we’ve built an inventory optimization model that:

Satisfies a 95% service level
Calculates the reorder point
Finds the order quantity that minimizes total cost

This is a powerful way to bring data-driven decisions to supply chain management!

Solving the Assignment Problem in Python

May 8, 2025

A Hands-On Example

The assignment problem is a classic optimization problem.
Suppose you have n workers and n jobs, and each worker has a different cost (or time, effort, etc.) associated with doing each job.
The goal is to assign each worker to exactly one job in a way that minimizes the total cost.

Mathematically, this is framed as minimizing the total cost:

$$
\min_{\sigma \in S_n} \sum_{i=1}^{n} C_{i, \sigma(i)}
$$

Where:

$C$ is the cost matrix.
$\sigma$ is a permutation of jobs assigned to workers.

Let’s walk through a concrete example and solve it using Python.

🧪 Example Scenario

We have 4 staff members and 4 jobs.
The cost matrix is as follows (rows = staff, columns = jobs):

	Job 1	Job 2	Job 3	Job 4
A	90	76	75	80
B	35	85	55	65
C	125	95	90	105
D	45	110	95	115

Our task is to assign each staff member to a job such that the total cost is minimized.

🔍 Step-by-Step: Python Code

We’ll use scipy.optimize.linear_sum_assignment, which implements the Hungarian Algorithm.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linear_sum_assignment

# Step 1: Define the cost matrix
cost_matrix = np.array([
    [90, 76, 75, 80],   # A
    [35, 85, 55, 65],   # B
    [125, 95, 90, 105], # C
    [45, 110, 95, 115]  # D
])

# Step 2: Solve the assignment problem
row_ind, col_ind = linear_sum_assignment(cost_matrix)
total_cost = cost_matrix[row_ind, col_ind].sum()

# St 3: Output the result
staff = ['A', 'B', 'C', 'D']
jobs = ['Job 1', 'Job 2', 'Job 3', 'Job 4']

print("Optimal Assignment:")
for i in range(len(row_ind)):
    print(f"  Staff {staff[row_ind[i]]} → {jobs[col_ind[i]]} (Cost: {cost_matrix[row_ind[i], col_ind[i]]})")

print(f"\nTotal Minimum Cost: {total_cost}")

💡 Code Explanation

Step 1: We define the cost_matrix, where each cell $(i, j)$ indicates the cost of assigning staff $i$ to job $j$.
Step 2: We use linear_sum_assignment() from SciPy to compute the optimal assignment.
- row_ind and col_ind give the optimal pairing.
Step 3: We print the assignments and compute the total cost.

📊 Visualization: Heatmap of Assignments

Let’s visualize the assignment on a heatmap where the assigned cells are highlighted.

# Step 4: Visualize with heatmap
plt.figure(figsize=(8, 6))
ax = sns.heatmap(cost_matrix, annot=True, fmt="d", cmap="YlGnBu", cbar=True,
                 xticklabels=jobs, yticklabels=staff)

# Highlight the assigned cells
for i in range(len(row_ind)):
    ax.add_patch(plt.Rectangle((col_ind[i], row_ind[i]), 1, 1, fill=False, edgecolor='red', lw=3))

plt.title(f"Optimal Assignment (Total Cost = {total_cost})")
plt.xlabel("Jobs")
plt.ylabel("Staff")
plt.tight_layout()
plt.show()

🧠 Result Interpretation

The red boxes in the heatmap show the optimal assignment.
Here’s what it might look like:

Optimal Assignment:
  Staff A → Job 4 (Cost: 80)
  Staff B → Job 3 (Cost: 55)
  Staff C → Job 2 (Cost: 95)
  Staff D → Job 1 (Cost: 45)

Total Minimum Cost: 275

This means if each staff member takes on the assigned job, the total cost will be minimized to 275.

📌 Conclusion

The assignment problem is a fundamental optimization challenge with applications in scheduling, resource allocation, and logistics.

Python’s scipy.optimize library makes it easy to solve real-world cases efficiently.

The heatmap visualization helps stakeholders intuitively understand the optimal decision.

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.