Optimizing Hospital Operating Room Schedules with Python and Linear Programming

April 12, 2025

Optimizing Operating Room Scheduling with Python: A Practical Approach

In healthcare systems, the efficient allocation of scarce medical resources like operating rooms (ORs) is crucial.
Poor scheduling can lead to underutilization of costly facilities or, worse, delays in critical surgeries.
In this blog, we’ll walk through a simplified example of OR scheduling using Python, model it as a Linear Programming (LP) problem, and visualize the results.

🏥 The Problem: Operating Room Scheduling

Imagine a hospital has $3$ operating rooms and needs to schedule $5$ surgeries.
Each surgery requires a different amount of time and has a different priority (urgency).
Our goal is to maximize total priority value while respecting the limited availability of operating rooms.

We assume:

Surgeries cannot be split.
Each operating room can handle at most one surgery.
Each surgery takes a fixed amount of time.
There are $8$ working hours per day per OR.

Let’s define the problem using the following notation:

Let $( x_{ij} )$ be a binary variable: $1$ if surgery $( i )$ is assigned to OR $( j )$, $0$ otherwise.
Each surgery $( i )$ has:
- Duration $( d_i )$ (in hours)
- Priority $( p_i )$
Total working hours per OR per day: $( H = 8 )$

📦 Data Setup and Model in Python

We’ll solve this using the PuLP library for Linear Programming.

!pip install pulp
import pulp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Define surgeries
surgeries = {
    'S1': {'duration': 3, 'priority': 10},
    'S2': {'duration': 4, 'priority': 8},
    'S3': {'duration': 2, 'priority': 6},
    'S4': {'duration': 5, 'priority': 7},
    'S5': {'duration': 1, 'priority': 5}
}

operating_rooms = ['OR1', 'OR2', 'OR3']
HOURS_PER_OR = 8

# Initialize LP problem
prob = pulp.LpProblem("OR_Scheduling", pulp.LpMaximize)

# Define decision variables
x = pulp.LpVariable.dicts("x", ((s, o) for s in surgeries for o in operating_rooms), cat="Binary")

# Objective: Maximize total priority
prob += pulp.lpSum(x[s, o] * surgeries[s]['priority'] for s in surgeries for o in operating_rooms)

# Constraint 1: Each surgery is assigned to at most one OR
for s in surgeries:
    prob += pulp.lpSum(x[s, o] for o in operating_rooms) <= 1

# Constraint 2: Each OR has a time budget of 8 hours
for o in operating_rooms:
    prob += pulp.lpSum(x[s, o] * surgeries[s]['duration'] for s in surgeries) <= HOURS_PER_OR

# Solve the problem
prob.solve()

✅ Explanation of the Code

Data Definition:
- We define $5$ surgeries with their durations and priority scores.
- There are $3$ available operating rooms.
Model Setup:
- We use pulp.LpProblem to define a maximization problem.
- Binary decision variables $( x_{ij} )$ are created to indicate assignments.
Objective Function:
- We maximize the total sum of priorities for scheduled surgeries.
Constraints:
- Each surgery can be assigned to at most one OR.
- Each OR can’t exceed its 8-hour daily limit.
Solving:
- We call .solve() to let PuLP find the optimal assignment.

📊 Visualizing the Results

Let’s show the assignments in a clear visual schedule.

# Parse results
schedule = []
for (s, o), var in x.items():
    if var.varValue == 1:
        schedule.append({
            'Surgery': s,
            'Operating Room': o,
            'Duration': surgeries[s]['duration'],
            'Priority': surgeries[s]['priority']
        })

df_schedule = pd.DataFrame(schedule)

# Add start and end time for visualization
df_schedule['Start'] = 0
df_schedule['End'] = df_schedule['Duration'].cumsum()

# Sort by OR
df_schedule.sort_values(by=['Operating Room', 'Priority'], ascending=[True, False], inplace=True)

# Plot Gantt chart
plt.figure(figsize=(10, 5))
sns.set_style("whitegrid")

for idx, row in df_schedule.iterrows():
    plt.barh(row['Operating Room'], row['Duration'], left=row['Start'], edgecolor='black', label=row['Surgery'])
    plt.text(row['Start'] + 0.5, row['Operating Room'], row['Surgery'], va='center', ha='left', color='white')

plt.xlabel("Hours")
plt.title("Operating Room Surgery Schedule")
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()

📈 Interpretation of the Output

The horizontal bar chart shows how surgeries are distributed across ORs.
You can see:
- Which surgeries got scheduled
- In which room
- For how long
The visual layout ensures no OR exceeds its 8-hour capacity.

🧠 Takeaways

Even a small scheduling problem can be effectively solved using Linear Programming.
PuLP in Python offers a clean way to formulate such optimization problems.
Visualization helps medical staff and administrators quickly interpret schedules.

In the real world, this problem could be expanded to consider:

Multiple days
Emergency surgery reservations
Surgeon availability
Pre-op/post-op dependencies

But even this simplified version illustrates how optimization can significantly improve healthcare delivery.

Optimizing Disaster Evacuation and Resource Allocation with Python

April 11, 2025

Here’s a explanation of a disaster evacuation and resource allocation optimization problem using Python.
This example is designed to run in Google Colaboratory.

🆘 Optimizing Evacuation Routes and Resource Allocation During Disasters with Python

In the face of natural disasters like earthquakes or floods, efficient evacuation and proper distribution of limited emergency resources are critical to saving lives.
In this blog post, we’ll walk through a practical example where we optimize:

Evacuation routes to shelters
Distribution of limited supplies (like food, water, and medicine)

We’ll solve this as a Linear Programming (LP) problem using Python’s PuLP library and visualize the results with matplotlib and networkx.

🧠 The Problem: Evacuation and Supply Allocation

Let’s consider a simplified city model with:

3 danger zones (where people need evacuation)
2 shelters (safe zones)
Limited roads connecting them
Each shelter has a capacity
A limited number of supplies that must be distributed optimally based on the number of evacuees

We’ll optimize for:

Minimum total travel distance
Feasible distribution of evacuees and resources

🧮 Mathematical Formulation

Let:

$( D = {d_1, d_2, d_3} )$ be danger zones
$( S = {s_1, s_2} )$ be shelters
$( x_{ij} )$: number of people evacuated from danger zone $( d_i )$ to shelter $( s_j )$
$( c_{ij} )$: cost (distance) from $( d_i )$ to $( s_j )$

Objective:

Minimize total cost:

$$
\min \sum_{i \in D} \sum_{j \in S} c_{ij} \cdot x_{ij}
$$

Subject to:

All people in each danger zone must be evacuated:
$$
\sum_{j \in S} x_{ij} = p_i \quad \forall i \in D
$$
Shelter capacity must not be exceeded:
$$
\sum_{i \in D} x_{ij} \leq C_j \quad \forall j \in S
$$

Where $( p_i )$ is the population in $( d_i )$ and $( C_j )$ is the capacity of $( s_j )$.

🧪 Let’s Code!

✅ Install PuLP

1	!pip install pulp

🧮 Python Code (Evacuation Optimization)

import pulp
import matplotlib.pyplot as plt
import networkx as nx

# Danger zones and shelter data
danger_zones = ['D1', 'D2', 'D3']
shelters = ['S1', 'S2']
population = {'D1': 100, 'D2': 150, 'D3': 80}
capacity = {'S1': 200, 'S2': 150}

# Distance cost matrix (from D to S)
costs = {
    ('D1', 'S1'): 10, ('D1', 'S2'): 20,
    ('D2', 'S1'): 30, ('D2', 'S2'): 10,
    ('D3', 'S1'): 20, ('D3', 'S2'): 30
}

# Define problem
prob = pulp.LpProblem("Evacuation_Optimization", pulp.LpMinimize)

# Decision variables
x = pulp.LpVariable.dicts("Evacuate", (danger_zones, shelters), lowBound=0, cat='Integer')

# Objective function
prob += pulp.lpSum([costs[(i, j)] * x[i][j] for i in danger_zones for j in shelters])

# Constraints: Each danger zone must evacuate all people
for i in danger_zones:
    prob += pulp.lpSum([x[i][j] for j in shelters]) == population[i]

# Constraints: Each shelter must not exceed its capacity
for j in shelters:
    prob += pulp.lpSum([x[i][j] for i in danger_zones]) <= capacity[j]

# Solve
prob.solve()

# Output results
for i in danger_zones:
    for j in shelters:
        print(f"Evacuate from {i} to {j}: {x[i][j].varValue} people")

print(f"Total evacuation cost: {pulp.value(prob.objective)}")

📊 Visualization

# Create graph for visualization
G = nx.DiGraph()

# Add nodes
G.add_nodes_from(danger_zones, bipartite=0)
G.add_nodes_from(shelters, bipartite=1)

# Add edges with flow values
edge_labels = {}
for i in danger_zones:
    for j in shelters:
        flow = x[i][j].varValue
        if flow > 0:
            G.add_edge(i, j, weight=flow)
            edge_labels[(i, j)] = f"{flow:.0f} ppl"

# Layout
pos = {
    'D1': (-1, 1), 'D2': (-1, 0), 'D3': (-1, -1),
    'S1': (1, 0.5), 'S2': (1, -0.5)
}

# Draw
plt.figure(figsize=(8, 5))
nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=2000, font_size=12)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.title("Evacuation Flow from Danger Zones to Shelters")
plt.show()

🔍 Results and Analysis

Evacuate from D1 to S1: 100.0 people
Evacuate from D1 to S2: 0.0 people
Evacuate from D2 to S1: 0.0 people
Evacuate from D2 to S2: 150.0 people
Evacuate from D3 to S1: 80.0 people
Evacuate from D3 to S2: 0.0 people
Total evacuation cost: 4100.0

The solver calculates the optimal evacuation plan minimizing total travel distance while respecting shelter capacities.
The graph clearly shows:

Which danger zones send evacuees to which shelters
The number of evacuees per route
How capacity constraints affect distribution

For example, if S1 is closer to D1 and D3, but its capacity is limited, the algorithm automatically diverts D2‘s population to S2 even if it’s farther, because that’s the optimal feasible plan.

🧯 Next Steps

Add resource distribution (food, medicine) using multi-objective optimization
Consider road congestion as dynamic costs
Extend to real geospatial data using libraries like geopandas

This example is just the tip of the iceberg.
Python and linear programming give us powerful tools to make life-saving decisions smarter and faster.

Optimizing Power Grid Design with Reliability-Aware Minimum Spanning Trees in Python

April 10, 2025

Let’s dive into a explanation of a reliable network design problem, solve it using Python, and visualize the result.
The example will focus on power grid design, where we aim to minimize cost while maintaining a high level of connectivity and reliability.
We’ll use NetworkX and Matplotlib for this.

🛠️ Designing a Reliable Power Grid with Python

Optimizing Cost and Reliability in a Network

In today’s world, designing infrastructure like power grids or communication networks requires not just cost-efficiency but also reliability.
What happens if a single node fails?
How do we ensure that power or data still flows smoothly?

This article tackles this question by modeling the Optimal Reliable Network Design problem.
We’ll explore a toy example and solve it with Python.

📘 The Problem Statement

Let’s assume we have a list of cities (nodes) that must be connected by power lines (edges).
Each connection has a cost and a reliability score (probability of functioning).
Our goal is to connect all cities (i.e., create a spanning network) while:

Minimizing total cost
Maximizing reliability, meaning if some connections fail, the network should still remain mostly connected.

We’ll use the K-Connected Minimum Spanning Tree concept — a spanning tree with redundancy.

🧮 Mathematical Formulation

Given:

A graph $( G = (V, E) )$
Cost $( c(e) )$ and reliability $( r(e) \in (0, 1] )$ for each edge $( e \in E )$

We want to find a subgraph $( T \subseteq G )$ such that:

$( T )$ is connected
The total cost $( \sum_{e \in T} c(e) )$ is minimized
The expected reliability is maximized:

$$
R(T) = \prod_{e \in T} r(e)
$$

We’ll simplify by selecting the Minimum Spanning Tree (MST) biased toward high reliability connections.

🐍 Python Code

Let’s implement this.

import networkx as nx
import matplotlib.pyplot as plt

# Create a graph
G = nx.Graph()

# Define nodes (cities)
nodes = ['A', 'B', 'C', 'D', 'E']
G.add_nodes_from(nodes)

# Define edges with cost and reliability
edges = [
    ('A', 'B', {'cost': 4, 'reliability': 0.95}),
    ('A', 'C', {'cost': 2, 'reliability': 0.90}),
    ('B', 'C', {'cost': 1, 'reliability': 0.99}),
    ('B', 'D', {'cost': 5, 'reliability': 0.85}),
    ('C', 'D', {'cost': 8, 'reliability': 0.80}),
    ('C', 'E', {'cost': 10, 'reliability': 0.70}),
    ('D', 'E', {'cost': 2, 'reliability': 0.95}),
]

G.add_edges_from([(u, v, attr) for u, v, attr in edges])

# Modify cost to incorporate reliability penalty
for u, v, attr in G.edges(data=True):
    reliability = attr['reliability']
    cost = attr['cost']
    # Adjusted cost: cost divided by reliability (to prefer high-reliability edges)
    attr['adjusted_cost'] = cost / reliability

# Build minimum spanning tree based on adjusted cost
mst = nx.minimum_spanning_tree(G, weight='adjusted_cost')

# Draw graph
pos = nx.spring_layout(G, seed=42)
plt.figure(figsize=(12, 6))

# Full graph
plt.subplot(121)
nx.draw(G, pos, with_labels=True, node_color='lightblue', edge_color='gray', node_size=800)
edge_labels = {(u, v): f"{d['cost']}, r={d['reliability']:.2f}" for u, v, d in G.edges(data=True)}
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.title("Original Graph (Cost, Reliability)")

# Minimum Spanning Tree
plt.subplot(122)
nx.draw(G, pos, with_labels=True, node_color='lightgreen', edge_color='lightgray', node_size=800)
nx.draw_networkx_edges(mst, pos, edge_color='red', width=2)
edge_labels_mst = {(u, v): f"{G[u][v]['adjusted_cost']:.2f}" for u, v in mst.edges()}
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels_mst)
plt.title("Minimum Spanning Tree (Adjusted Cost)")

plt.tight_layout()
plt.show()

# Calculate total cost and reliability
total_cost = sum(G[u][v]['cost'] for u, v in mst.edges())
total_reliability = 1
for u, v in mst.edges():
    total_reliability *= G[u][v]['reliability']

print(f"Total MST Cost: {total_cost}")
print(f"Expected Network Reliability: {total_reliability:.4f}")

🔍 Code Breakdown

Graph Construction: We create a graph where each edge has cost and reliability.
Adjusted Cost: We penalize unreliable connections by dividing cost by reliability.
This favors low-cost, high-reliability edges.
MST Generation: We use NetworkX’s built-in minimum_spanning_tree() function with our adjusted cost.
Visualization: Two side-by-side graphs — one for the original network, and one showing the selected MST.
Reliability Estimation: We assume independence of edge failures and multiply the probabilities.

📊 Results & Interpretation

Total MST Cost: 10
Expected Network Reliability: 0.7195

The algorithm produces a network that:

Connects all cities
Minimizes the cost while penalizing unreliable connections
Shows how adjusting the weights influences the structure of the network

The MST may avoid the lowest-cost edge if it has poor reliability.
This showcases how simple weighting schemes can enforce more robust designs.

✅ Conclusion

This small example shows how multi-objective network optimization can be approached using Python.
In real-world problems, constraints like minimum node degree, backup paths, or capacity limits can be incorporated.

Optimizing Investment Portfolios

April 9, 2025

Solving the Capital Allocation Problem with Python

Today, I’m excited to explore an interesting financial mathematics problem: the capital allocation problem.
This problem is fundamental in portfolio management, where we need to decide how to distribute funds across different investment opportunities to achieve optimal returns while managing risk.

Let’s dive into a concrete example and solve it using Python.
I’ll walk you through the entire process, from mathematical formulation to visualization of results.

Problem Statement

Imagine we have $10,000 to invest in 4 different assets: stocks, bonds, real estate, and commodities.
Each asset has its own expected return, volatility (risk), and correlation with other assets.
Our goal is to find the optimal allocation that maximizes the Sharpe ratio (return per unit of risk).

Mathematical Formulation

The Sharpe ratio is defined as:

$$S = \frac{E[R] - R_f}{\sigma}$$

Where:

$E[R]$ is the expected return of the portfolio
$R_f$ is the risk-free rate
$\sigma$ is the standard deviation (volatility) of the portfolio

For a portfolio with weights $w_i$ for each asset $i$, the expected return is:

$$E[R_p] = \sum_{i=1}^{n} w_i E[R_i]$$

And the portfolio variance is:

$$\sigma_p^2 = \sum_{i=1}^{n} \sum_{j=1}^{n} w_i w_j \sigma_i \sigma_j \rho_{ij}$$

Where $\rho_{ij}$ is the correlation between assets $i$ and $j$.

Let’s implement this in Python:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
from matplotlib.cm import get_cmap

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

# Define our assets
assets = ['Stocks', 'Bonds', 'Real Estate', 'Commodities']
n_assets = len(assets)

# Expected returns (annual)
expected_returns = np.array([0.10, 0.05, 0.08, 0.12])

# Volatilities (annual)
volatilities = np.array([0.20, 0.08, 0.15, 0.25])

# Create a correlation matrix
# This is a symmetric matrix with 1s on the diagonal
correlation_matrix = np.array([
    [1.00, 0.20, 0.50, 0.30],
    [0.20, 1.00, 0.30, 0.10],
    [0.50, 0.30, 1.00, 0.40],
    [0.30, 0.10, 0.40, 1.00]
])

# Calculate the covariance matrix
covariance_matrix = np.zeros((n_assets, n_assets))
for i in range(n_assets):
    for j in range(n_assets):
        covariance_matrix[i, j] = correlation_matrix[i, j] * volatilities[i] * volatilities[j]

# Visualize the correlation matrix
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', xticklabels=assets, yticklabels=assets)
plt.title('Asset Correlation Matrix')
plt.tight_layout()
plt.show()

# Risk-free rate
risk_free_rate = 0.02

# Function to calculate portfolio return and volatility
def portfolio_performance(weights):
    # Returns
    returns = np.sum(weights * expected_returns)
    
    # Volatility
    volatility = np.sqrt(np.dot(weights.T, np.dot(covariance_matrix, weights)))
    
    return returns, volatility

# Function to minimize (negative Sharpe ratio)
def negative_sharpe_ratio(weights):
    returns, volatility = portfolio_performance(weights)
    sharpe_ratio = (returns - risk_free_rate) / volatility
    return -sharpe_ratio

# Constraints: weights sum to 1
constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - 1}

# Bounds: weights between 0 and 1 (no short selling)
bounds = tuple((0, 1) for _ in range(n_assets))

# Initial guess: equal weights
initial_weights = np.array([1/n_assets] * n_assets)

# Optimize
result = minimize(negative_sharpe_ratio, initial_weights, method='SLSQP', 
                 bounds=bounds, constraints=constraints)

# Extract the optimal weights
optimal_weights = result['x']

# Calculate the performance of the optimal portfolio
optimal_returns, optimal_volatility = portfolio_performance(optimal_weights)
optimal_sharpe = (optimal_returns - risk_free_rate) / optimal_volatility

print("Optimal Portfolio Weights:")
for asset, weight in zip(assets, optimal_weights):
    print(f"{asset}: {weight:.4f} ({weight*100:.2f}%)")
    
print(f"\nOptimal Portfolio Performance:")
print(f"Expected Annual Return: {optimal_returns:.4f} ({optimal_returns*100:.2f}%)")
print(f"Annual Volatility: {optimal_volatility:.4f} ({optimal_volatility*100:.2f}%)")
print(f"Sharpe Ratio: {optimal_sharpe:.4f}")

# Calculate total investment amounts (assuming $10,000 total)
total_investment = 10000
investment_amounts = optimal_weights * total_investment

print(f"\nCapital Allocation for ${total_investment}:")
for asset, amount in zip(assets, investment_amounts):
    print(f"{asset}: ${amount:.2f}")

# Monte Carlo simulation to generate random portfolios
num_portfolios = 10000
results = np.zeros((3, num_portfolios))
all_weights = np.zeros((num_portfolios, n_assets))

for i in range(num_portfolios):
    # Generate random weights
    weights = np.random.random(n_assets)
    weights = weights / np.sum(weights)
    all_weights[i, :] = weights
    
    # Calculate portfolio return and volatility
    portfolio_return, portfolio_volatility = portfolio_performance(weights)
    
    # Store results
    results[0, i] = portfolio_return
    results[1, i] = portfolio_volatility
    results[2, i] = (portfolio_return - risk_free_rate) / portfolio_volatility

# Plot the efficient frontier
plt.figure(figsize=(12, 8))
plt.scatter(results[1, :], results[0, :], c=results[2, :], cmap='viridis', 
            s=10, alpha=0.3, marker='o')

# Plot the optimal portfolio
plt.scatter(optimal_volatility, optimal_returns, c='red', s=100, marker='*', 
            label=f'Optimal Portfolio (Sharpe: {optimal_sharpe:.4f})')

# Add labels and legend
plt.colorbar(label='Sharpe Ratio')
plt.xlabel('Volatility (Standard Deviation)')
plt.ylabel('Expected Return')
plt.title('Portfolio Optimization: Risk vs Return')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Plot the optimal portfolio allocation
plt.figure(figsize=(12, 6))
colors = get_cmap('tab10').colors[:n_assets]

# Bar chart
plt.subplot(1, 2, 1)
plt.bar(assets, optimal_weights, color=colors)
plt.title('Optimal Portfolio Weights')
plt.ylabel('Weight')
plt.xticks(rotation=45)

# Pie chart
plt.subplot(1, 2, 2)
plt.pie(optimal_weights, labels=assets, autopct='%1.1f%%', colors=colors)
plt.title('Optimal Portfolio Allocation')

plt.tight_layout()
plt.show()

# Test some alternative allocations
alternative_allocations = [
    {'name': 'Equal Weight', 'weights': np.array([0.25, 0.25, 0.25, 0.25])},
    {'name': 'Stocks Heavy', 'weights': np.array([0.70, 0.10, 0.10, 0.10])},
    {'name': 'Bonds Heavy', 'weights': np.array([0.10, 0.70, 0.10, 0.10])},
    {'name': 'Optimal', 'weights': optimal_weights}
]

# Calculate performance metrics for each allocation
performance_data = []
for alloc in alternative_allocations:
    returns, volatility = portfolio_performance(alloc['weights'])
    sharpe = (returns - risk_free_rate) / volatility
    performance_data.append({
        'Allocation': alloc['name'],
        'Return': returns * 100,  # Convert to percentage
        'Volatility': volatility * 100,  # Convert to percentage
        'Sharpe Ratio': sharpe
    })

# Convert to DataFrame for easier visualization
performance_df = pd.DataFrame(performance_data)

# Create a comparison chart
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# Return comparison
axes[0].bar(performance_df['Allocation'], performance_df['Return'], color='skyblue')
axes[0].set_title('Expected Annual Return (%)')
axes[0].set_ylabel('Return (%)')
axes[0].set_ylim(0, max(performance_df['Return']) * 1.2)
axes[0].grid(axis='y', alpha=0.3)
axes[0].set_xticklabels(performance_df['Allocation'], rotation=45)

# Volatility comparison
axes[1].bar(performance_df['Allocation'], performance_df['Volatility'], color='salmon')
axes[1].set_title('Annual Volatility (%)')
axes[1].set_ylabel('Volatility (%)')
axes[1].set_ylim(0, max(performance_df['Volatility']) * 1.2)
axes[1].grid(axis='y', alpha=0.3)
axes[1].set_xticklabels(performance_df['Allocation'], rotation=45)

# Sharpe Ratio comparison
axes[2].bar(performance_df['Allocation'], performance_df['Sharpe Ratio'], color='lightgreen')
axes[2].set_title('Sharpe Ratio')
axes[2].set_ylabel('Sharpe Ratio')
axes[2].set_ylim(0, max(performance_df['Sharpe Ratio']) * 1.2)
axes[2].grid(axis='y', alpha=0.3)
axes[2].set_xticklabels(performance_df['Allocation'], rotation=45)

plt.tight_layout()
plt.show()

Code Explanation

Let’s break down what the code does:

Setup: We first define our assets, expected returns, volatilities, and correlation matrix.
These are the inputs to our optimization problem.
Covariance Matrix: We calculate the covariance matrix from the correlation matrix and volatilities.
This is essential for calculating portfolio risk.
Portfolio Performance: We create functions to calculate the expected return and volatility of a portfolio given certain weights.
Optimization: We use scipy.optimize.minimize to find the weights that maximize the Sharpe ratio (or minimize the negative Sharpe ratio).
We set constraints (weights sum to 1) and bounds (no short selling, so weights are between $0$ and $1$).
Monte Carlo Simulation: We generate $10,000$ random portfolios to visualize the efficient frontier.
Performance Comparison: We compare our optimal portfolio with some alternative allocations (equal weights, stocks-heavy, bonds-heavy).

Results and Visualization

The optimization process identifies the optimal portfolio weights that maximize the Sharpe ratio.
Let’s look at the key findings:

Optimal Allocation: The optimization suggests a diversified portfolio with specific weights for each asset class.

Optimal Portfolio Weights:
Stocks: 0.1700 (17.00%)
Bonds: 0.5481 (54.81%)
Real Estate: 0.1113 (11.13%)
Commodities: 0.1706 (17.06%)

Optimal Portfolio Performance:
Expected Annual Return: 0.0738 (7.38%)
Annual Volatility: 0.0927 (9.27%)
Sharpe Ratio: 0.5802

Capital Allocation for $10000:
Stocks: $1700.10
Bonds: $5480.85
Real Estate: $1113.21
Commodities: $1705.84

Risk-Return Profile: The efficient frontier plot shows the relationship between risk and return for different portfolios. The optimal portfolio (marked with a red star) provides the best risk-adjusted return.

Capital Allocation: For our $10,000 investment, we can see exactly how much to allocate to each asset class.

Comparison: The bar charts compare the optimal allocation with alternative strategies in terms of expected return, volatility, and Sharpe ratio.

Insights and Conclusions

This example demonstrates the power of modern portfolio theory and optimization techniques in solving capital allocation problems.
The key insights are:

Diversification Benefits: The optimal portfolio is well-diversified across all asset classes, taking advantage of the correlation structure to reduce overall risk.
Risk-Return Tradeoff: There’s a clear tradeoff between risk and return, but optimization helps us find the sweet spot that maximizes risk-adjusted returns.
Practical Application: With this approach, investors can make data-driven decisions about capital allocation rather than relying on intuition alone.

In practice, this approach would be enhanced with more sophisticated models, historical data analysis, and forward-looking scenarios.
However, the fundamental principles demonstrated here remain the same across more complex implementations.

The capital allocation problem is a perfect example of how mathematical optimization can be applied to real-world financial decisions, helping individuals and institutions make better investment choices.

Optimizing Machine-to-Task Assignments in Python with Visualization

April 8, 2025

🛠️ Solving the Machine Assignment Problem with Python: A Step-by-Step Guide

In the world of operations research and optimization, the machine assignment problem is a classic example of how to allocate resources (machines) to tasks (jobs) in the most efficient way possible.
The goal is to minimize the total cost or maximize the total profit associated with the assignments.

Let’s walk through a concrete example and solve it using Python on Google Colab.
We’ll use scipy.optimize to get an optimal solution, and matplotlib to visualize the result.

📘 Problem Statement

Suppose we have 3 machines and 3 tasks.
Each machine can handle only one task, and each task must be assigned to exactly one machine.
The cost matrix below shows the cost of assigning machine $( i )$ to task $( j )$:

$$
\text{Cost Matrix} =
\begin{bmatrix}
9 & 2 & 7 \
6 & 4 & 3 \
5 & 8 & 1 \
\end{bmatrix}
$$

Our goal is to find the optimal assignment of machines to tasks that minimizes the total cost.

🔍 Mathematical Formulation

This is a Linear Sum Assignment Problem, also known as the Assignment Problem.

Let:

$( C_{ij} )$: cost of assigning machine $( i )$ to task $( j )$
$( x_{ij} )$: binary variable, $1$ if machine $( i )$ is assigned to task $( j )$, $0$ otherwise

The objective function:

$$
\min \sum_{i=1}^{n} \sum_{j=1}^{n} C_{ij} \cdot x_{ij}
$$

Subject to:

$$
\sum_{j=1}^{n} x_{ij} = 1 \quad \forall i \quad \text{(each machine assigned to one task)}
$$
$$
\sum_{i=1}^{n} x_{ij} = 1 \quad \forall j \quad \text{(each task assigned to one machine)}
$$

🧠 Python Implementation

We will use scipy.optimize.linear_sum_assignment to solve the problem.

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

# Define the cost matrix
cost_matrix = np.array([
    [9, 2, 7],
    [6, 4, 3],
    [5, 8, 1]
])

# Solve the assignment problem
row_ind, col_ind = linear_sum_assignment(cost_matrix)

# Output the assignments and total cost
total_cost = cost_matrix[row_ind, col_ind].sum()

# Print results
print("Optimal assignment (Machine -> Task):")
for r, c in zip(row_ind, col_ind):
    print(f"Machine {r} -> Task {c} (Cost: {cost_matrix[r][c]})")
print(f"\nTotal Minimum Cost: {total_cost}")

🧩 Code Explanation

cost_matrix: A 3x3 matrix representing the cost for each machine-task assignment.
linear_sum_assignment(): Uses the Hungarian algorithm to find the optimal assignment.
row_ind, col_ind: Arrays that indicate which machine is assigned to which task.
We then compute the total cost from the selected assignments.

📊 Visualizing the Assignment

Let’s make a heatmap to visualize the assignments and costs.

import seaborn as sns

# Create a heatmap to show costs
plt.figure(figsize=(6, 5))
sns.heatmap(cost_matrix, annot=True, fmt="d", cmap="YlGnBu", cbar=False)

# Highlight the chosen assignments
for i, j in zip(row_ind, col_ind):
    plt.text(j + 0.2, i + 0.4, '✓', fontsize=18, color='red')

plt.title("Machine Assignment Cost Matrix")
plt.xlabel("Tasks")
plt.ylabel("Machines")
plt.xticks(np.arange(3) + 0.5, [f"Task {i}" for i in range(3)])
plt.yticks(np.arange(3) + 0.5, [f"Machine {i}" for i in range(3)])
plt.grid(False)
plt.show()

📌 Result Interpretation

The output of our code would look like this:

Optimal assignment (Machine -> Task):
Machine 0 -> Task 1 (Cost: 2)
Machine 1 -> Task 0 (Cost: 6)
Machine 2 -> Task 2 (Cost: 1)

Total Minimum Cost: 9

In the heatmap, you’ll see red check marks ✓ on the optimal machine-task pairs.

✅ Conclusion

We’ve successfully tackled the machine assignment problem using Python.
This type of problem is highly applicable in manufacturing, logistics, scheduling, and project management.
The scipy.optimize library offers a powerful and easy-to-use tool for solving these types of optimization challenges efficiently.

Next time you need to assign resources optimally — whether they’re machines, people, or delivery trucks — remember this approach!

Optimizing Delivery Routes：Solving Vehicle Routing Problems with Python

April 7, 2025

Today, I want to share a practical approach to tackling the Vehicle Routing Problem (VRP) - one of the most fascinating challenges in logistics and operations research. If you’ve ever wondered how delivery companies optimize their routes to serve multiple customers with limited vehicles, you’re about to find out!

Understanding the Vehicle Routing Problem

The Vehicle Routing Problem involves finding optimal routes for a fleet of vehicles to deliver goods to a set of customers. The objective is typically to minimize the total distance traveled while satisfying various constraints like vehicle capacity, time windows, etc.

Mathematically, the basic VRP can be expressed as:

$$\min \sum_{i=0}^{n} \sum_{j=0}^{n} \sum_{k=1}^{m} c_{ij} x_{ijk}$$

Subject to constraints like:
$$\sum_{k=1}^{m} \sum_{j=0}^{n} x_{ijk} = 1 \quad \forall i \in {1, 2, …, n}$$
$$\sum_{j=1}^{n} q_j \sum_{i=0}^{n} x_{ijk} \leq Q_k \quad \forall k \in {1, 2, …, m}$$

Where:

$c_{ij}$ is the cost (or distance) from location $i$ to $j$
$x_{ijk}$ equals 1 if vehicle $k$ travels from $i$ to $j$, and 0 otherwise
$q_j$ is the demand at customer $j$
$Q_k$ is the capacity of vehicle $k$

Let’s implement a solution to a concrete example!

A Practical VRP Example

Consider a company that needs to deliver packages to 10 customers from a central depot. The company has 3 vehicles, each with a capacity of 10 units. Each customer has a specific demand, and we want to minimize the total distance traveled.

import numpy as np
import matplotlib.pyplot as plt
import random
from itertools import permutations
from math import sqrt

class VRPSolver:
    def __init__(self, depot, customers, demands, vehicle_capacity, num_vehicles):
        """
        Initialize the VRP solver with problem parameters
        
        Parameters:
        depot: tuple (x, y) - coordinates of the depot
        customers: list of tuples [(x1, y1), (x2, y2), ...] - coordinates of customers
        demands: list [d1, d2, ...] - demand of each customer
        vehicle_capacity: int - capacity of each vehicle
        num_vehicles: int - number of available vehicles
        """
        self.depot = depot
        self.customers = customers
        self.demands = demands
        self.vehicle_capacity = vehicle_capacity
        self.num_vehicles = num_vehicles
        self.distance_matrix = self._compute_distance_matrix()
        
    def _compute_distance_matrix(self):
        """Compute the distance matrix between all locations (depot + customers)"""
        all_locations = [self.depot] + self.customers
        n = len(all_locations)
        dist_matrix = np.zeros((n, n))
        
        for i in range(n):
            for j in range(n):
                if i != j:
                    x1, y1 = all_locations[i]
                    x2, y2 = all_locations[j]
                    dist_matrix[i, j] = sqrt((x2 - x1)**2 + (y2 - y1)**2)
                    
        return dist_matrix
    
    def _calculate_route_cost(self, route):
        """Calculate the total distance of a route"""
        total_distance = 0
        prev_location = 0  # Start at depot (index 0)
        
        for customer in route:
            total_distance += self.distance_matrix[prev_location, customer]
            prev_location = customer
            
        # Return to depot
        total_distance += self.distance_matrix[prev_location, 0]
        
        return total_distance
    
    def _is_route_feasible(self, route):
        """Check if a route satisfies capacity constraints"""
        total_demand = sum(self.demands[customer-1] for customer in route)
        return total_demand <= self.vehicle_capacity
    
    def solve_greedy(self):
        """Solve VRP using a greedy algorithm"""
        unassigned = list(range(1, len(self.customers) + 1))
        routes = [[] for _ in range(self.num_vehicles)]
        
        # Assign customers to vehicles until all are assigned
        for vehicle_idx in range(self.num_vehicles):
            current_location = 0  # Depot
            remaining_capacity = self.vehicle_capacity
            
            while unassigned:
                best_next = None
                best_distance = float('inf')
                
                # Find the nearest unassigned customer that fits in the vehicle
                for customer in unassigned:
                    if self.demands[customer-1] <= remaining_capacity:
                        distance = self.distance_matrix[current_location, customer]
                        if distance < best_distance:
                            best_distance = distance
                            best_next = customer
                
                if best_next is None:
                    # No more customers can fit in this vehicle
                    break
                    
                # Assign the customer to this vehicle
                routes[vehicle_idx].append(best_next)
                remaining_capacity -= self.demands[best_next-1]
                current_location = best_next
                unassigned.remove(best_next)
                
            if not unassigned:
                break
        
        return routes
    
    def solve_heuristic(self):
        """Solve VRP using a simple heuristic with local search"""
        # Start with a greedy solution
        routes = self.solve_greedy()
        
        # Try to improve each route using 2-opt local search
        for i in range(len(routes)):
            if len(routes[i]) >= 2:
                routes[i] = self._improve_route_2opt(routes[i])
        
        return routes
    
    def _improve_route_2opt(self, route):
        """Improve a route using 2-opt local search"""
        improved = True
        best_route = route.copy()
        best_cost = self._calculate_route_cost(best_route)
        
        while improved:
            improved = False
            
            for i in range(len(route)):
                for j in range(i + 2, len(route)):
                    # Skip if indices would make an invalid swap
                    if j - i == 1:
                        continue
                        
                    # Create a new route with a 2-opt swap
                    new_route = route[:i + 1] + route[i + 1:j + 1][::-1] + route[j + 1:]
                    new_cost = self._calculate_route_cost(new_route)
                    
                    if new_cost < best_cost:
                        best_route = new_route
                        best_cost = new_cost
                        improved = True
                        
            route = best_route.copy()
        
        return best_route
    
    def visualize_solution(self, routes):
        """Visualize the VRP solution"""
        plt.figure(figsize=(10, 8))
        
        # Plot depot
        plt.scatter([self.depot[0]], [self.depot[1]], c='red', s=100, marker='s', label='Depot')
        
        # Plot customers
        customer_x = [c[0] for c in self.customers]
        customer_y = [c[1] for c in self.customers]
        plt.scatter(customer_x, customer_y, c='blue', s=50, label='Customers')
        
        # Add customer numbers
        for i, (x, y) in enumerate(self.customers):
            plt.annotate(f"{i+1}\n(d:{self.demands[i]})", (x, y), xytext=(5, 5), 
                         textcoords='offset points')
        
        # Plot routes with different colors
        colors = ['green', 'purple', 'orange', 'brown', 'pink', 'gray', 'olive', 'cyan']
        for i, route in enumerate(routes):
            if not route:
                continue
                
            # Get coordinates for the route (including depot at start and end)
            route_x = [self.depot[0]]
            route_y = [self.depot[1]]
            
            for customer in route:
                route_x.append(self.customers[customer-1][0])
                route_y.append(self.customers[customer-1][1])
                
            # Return to depot
            route_x.append(self.depot[0])
            route_y.append(self.depot[1])
            
            plt.plot(route_x, route_y, c=colors[i % len(colors)], label=f'Vehicle {i+1}')
        
        plt.title("Vehicle Routing Problem Solution")
        plt.xlabel("X Coordinate")
        plt.ylabel("Y Coordinate")
        plt.legend()
        plt.grid(True)
        plt.axis('equal')
        plt.tight_layout()
        
        return plt

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

# Problem setup
depot = (0, 0)
num_customers = 10
customer_locations = [(random.uniform(-10, 10), random.uniform(-10, 10)) for _ in range(num_customers)]
customer_demands = [random.randint(1, 5) for _ in range(num_customers)]
vehicle_capacity = 10
num_vehicles = 3

# Create and solve the VRP
vrp = VRPSolver(depot, customer_locations, customer_demands, vehicle_capacity, num_vehicles)
solution = vrp.solve_heuristic()

# Print the solution
print("Solution:")
total_distance = 0
for i, route in enumerate(solution):
    route_distance = vrp._calculate_route_cost(route)
    total_distance += route_distance
    route_demand = sum(vrp.demands[customer-1] for customer in route)
    print(f"Vehicle {i+1}: {route} (Distance: {route_distance:.2f}, Load: {route_demand}/{vehicle_capacity})")
print(f"Total distance: {total_distance:.2f}")

# Visualize the solution
plt = vrp.visualize_solution(solution)
plt.show()

Code Explanation

Let’s walk through the key components of this solution:

1. `VRPSolver` Class

The heart of our solution is the VRPSolver class, which encapsulates all the logic for solving the Vehicle Routing Problem:

Initialization: We set up the problem parameters (depot location, customer locations, demands, vehicle capacity, and number of vehicles) and compute a distance matrix.
Distance Calculation: The _compute_distance_matrix() method calculates the Euclidean distance between all locations (depot and customers).
Route Evaluation: _calculate_route_cost() computes the total distance of a route, and _is_route_feasible() checks if a route satisfies capacity constraints.

2. Solution Approaches

The code implements two solution methods:

Greedy Algorithm (solve_greedy()): This approach assigns customers to vehicles one by one, always choosing the nearest unassigned customer that fits within the vehicle’s remaining capacity.
Local Search Heuristic (solve_heuristic()): This method starts with the greedy solution and then applies a 2-opt local search to improve each route. The 2-opt algorithm works by swapping two edges in a route to see if a shorter path can be found.

3. Visualization

The visualize_solution() method creates a clear representation of the solution, showing:

The depot (red square)
Customer locations (blue circles) with their ID and demand
Vehicle routes with different colors

Results Analysis

When we run the code, we get both a textual and visual representation of the solution:

Solution:
Vehicle 1: [8, 3, 9, 6, 7] (Distance: 36.89, Load: 10/10)
Vehicle 2: [2, 5, 1] (Distance: 26.24, Load: 8/10)
Vehicle 3: [4, 10] (Distance: 25.35, Load: 7/10)
Total distance: 88.48

The solution divides the 10 customers among 3 vehicles, ensuring that each vehicle doesn’t exceed its capacity of 10 units. Let’s analyze what makes this a good solution:

Balanced Load: Each vehicle is using a significant portion of its capacity (one vehicles at 100%, one vehicles at 80% and one at 70%).
Route Efficiency: The routes minimize unnecessary travel by grouping customers that are geographically close.
Total Distance: The algorithm has found a solution with a total distance of 88.48 units.

The visualization shows us how the vehicles are routed. Each colored line represents a different vehicle’s route, starting from the depot, visiting assigned customers, and returning to the depot.

Optimization Opportunities

Our solution uses relatively simple heuristics. In a real-world scenario, you might consider:

More Advanced Algorithms: Such as simulated annealing, genetic algorithms, or exact methods like branch-and-cut for smaller instances.
Additional Constraints: Like time windows, multiple depots, or heterogeneous vehicles.
Dynamic Routing: Handling real-time changes, such as new orders or traffic conditions.

Conclusion

Vehicle Routing Problems are at the core of modern logistics operations. Even with our relatively simple implementation, we can see how mathematical modeling and algorithms can help optimize delivery routes, saving costs and time.

The next time you receive a package delivery, think about the complex optimization problem that was solved behind the scenes to get that package to your doorstep efficiently!

Designing Cost-Effective Networks with Python：A Visual Guide to Minimum Spanning Trees

April 6, 2025

Let’s walk through a real-world example of a network design problem, solve it using Python, and visualize the results.

🧠 Solving a Network Design Problem with Python

Network design problems are everywhere — from planning internet infrastructure to organizing warehouse logistics.
Today, we’ll solve a simplified version using Python, and visualize our solution to really bring it to life.

📌 Problem Statement: Minimum Cost Network Design

Imagine you’re building a network to connect a number of office buildings in a city.
You want to connect all the buildings using network cables.
Your goal is:

Connect all buildings with the minimum total cost, where the cost is based on the distance between them.

This is a classic Minimum Spanning Tree (MST) problem in graph theory.
We’ll use Prim’s Algorithm to solve it.

📊 Sample Data

Let’s say we have 6 buildings.
The costs (distances) between them are represented in the following adjacency matrix:

	A	B	C	D	E	F
A	0	3	0	0	6	5
B	3	0	1	0	0	4
C	0	1	0	6	0	4
D	0	0	6	0	8	5
E	6	0	0	8	0	2
F	5	4	4	5	2	0

🔢 Mathematical Formulation

We want to find a subgraph $( T = (V, E_T) )$ of the graph $( G = (V, E) )$ such that:

$( T )$ is connected,
$( T )$ spans all nodes (buildings),
$( \sum_{e \in E_T} w(e) )$ is minimized (total cost is minimum).

🐍 Python Code: Prim’s Algorithm

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

# Define the adjacency matrix
cost_matrix = np.array([
    [0, 3, 0, 0, 6, 5],
    [3, 0, 1, 0, 0, 4],
    [0, 1, 0, 6, 0, 4],
    [0, 0, 6, 0, 8, 5],
    [6, 0, 0, 8, 0, 2],
    [5, 4, 4, 5, 2, 0]
])

labels = ['A', 'B', 'C', 'D', 'E', 'F']
num_nodes = len(cost_matrix)

# Prim's Algorithm
def prim_mst(cost_matrix):
    selected_nodes = [False] * num_nodes
    selected_nodes[0] = True
    edges = []

    for _ in range(num_nodes - 1):
        min_cost = float('inf')
        a = b = -1
        for i in range(num_nodes):
            if selected_nodes[i]:
                for j in range(num_nodes):
                    if not selected_nodes[j] and cost_matrix[i][j]:
                        if cost_matrix[i][j] < min_cost:
                            min_cost = cost_matrix[i][j]
                            a, b = i, j
        if a != -1 and b != -1:
            edges.append((a, b, cost_matrix[a][b]))
            selected_nodes[b] = True

    return edges

# Compute MST
mst_edges = prim_mst(cost_matrix)

🧠 Code Explanation

We define a cost (distance) matrix where 0 means “no direct connection.”
We start from node 0 (building A), and at each step, add the minimum cost edge connecting a new node.
The result is a list of edges forming the MST.

📈 Visualizing the Network

# Create graph and MST subgraph
G = nx.Graph()
MST = nx.Graph()

# Add all edges for base graph
for i in range(num_nodes):
    for j in range(i + 1, num_nodes):
        if cost_matrix[i][j] != 0:
            G.add_edge(labels[i], labels[j], weight=cost_matrix[i][j])

# Add only MST edges
for i, j, w in mst_edges:
    MST.add_edge(labels[i], labels[j], weight=w)

# Position nodes in a circular layout
pos = nx.circular_layout(G)

# Plot full graph with light edges
plt.figure(figsize=(10, 6))
nx.draw_networkx(G, pos, with_labels=True, node_color='lightblue', edge_color='gray', node_size=1200)
nx.draw_networkx_edges(MST, pos, edge_color='red', width=3)
edge_labels = nx.get_edge_attributes(MST, 'weight')
nx.draw_networkx_edge_labels(MST, pos, edge_labels=edge_labels)

plt.title("Minimum Spanning Tree (in Red) over Full Network")
plt.axis('off')
plt.show()

📌 Final Result

The red edges in the graph represent the optimal set of connections (the MST).
These are the links you’d install in your network to connect all buildings with the least amount of cabling cost.

The MST edges and costs are:

1 2	for u, v, w in mst_edges: print(f"{labels[u]} - {labels[v]} : {w}")

Output:

A - B : 3  
B - C : 1  
B - F : 4  
F - E : 2  
F - D : 5

✅ Total cost: $(3 + 1 + 4 + 2 + 5 = \boxed{15})$

🧩 Conclusion

This example shows how Python can be a powerful tool for network design optimization, even in introductory scenarios.
We used:

Graph theory (MST)
A real algorithm (Prim’s)
Visualization (NetworkX + Matplotlib)

Optimizing Material Cuts：Solving the Cutting Stock Problem with Python

April 5, 2025

Solving the Cutting Stock Problem with Python 🧠✂️

The Cutting Stock Problem (CSP) is a classic optimization issue that appears in many industries, especially manufacturing.
The challenge is to minimize waste when cutting raw materials into smaller pieces of specified lengths.

Let’s walk through a concrete example and solve it using Python.
We’ll also visualize the results to make everything crystal clear.

🧩 Problem Statement

Imagine a company that cuts steel rods of length $100$ units into smaller pieces to fulfill customer orders.

Here’s the customer demand:

Length Needed	Quantity
20	48
45	35
50	24

The company wants to minimize the number of $100$-unit rods used while fulfilling this demand.
This is a classic integer linear programming problem, often solved using column generation techniques.
But we’ll approach it in a simplified way using PuLP, a Python library for linear programming.

📌 Mathematical Formulation

Let’s define:

Raw stock length: $( L = 100 )$
Demand: $( d_i ) units of size ( l_i )$

We generate possible cutting patterns that fit within $100$ units.
Each pattern becomes a variable in the optimization model.

Let:

$( x_j )$: number of times pattern $( j )$ is used
$( a_{ij} )$: number of pieces of type $( i )$ in pattern $( j )$

We aim to:

$$
\min \sum_j x_j
$$

Subject to:

$$
\sum_j a_{ij} x_j \geq d_i \quad \text{for all } i
$$

$$
x_j \in \mathbb{Z}_{\geq 0}
$$

🐍 Python Implementation

Let’s jump into the code.
We’ll use PuLP to solve this as an integer program and matplotlib to visualize the results.

# Install PuLP if needed
!pip install pulp

import pulp
import itertools
import matplotlib.pyplot as plt

# Parameters
stock_length = 100
demand = {20: 48, 45: 35, 50: 24}
sizes = list(demand.keys())

# Generate all valid cutting patterns
def generate_patterns(sizes, stock_length):
    patterns = []
    for counts in itertools.product(*(range(stock_length // s + 1) for s in sizes)):
        total = sum(s * c for s, c in zip(sizes, counts))
        if 0 < total <= stock_length:
            patterns.append(counts)
    return patterns

patterns = generate_patterns(sizes, stock_length)

# Define the optimization problem
prob = pulp.LpProblem("CuttingStockProblem", pulp.LpMinimize)

# Define variables: how many times to use each pattern
x = [pulp.LpVariable(f"x_{i}", lowBound=0, cat='Integer') for i in range(len(patterns))]

# Objective: minimize number of stock rods used
prob += pulp.lpSum(x)

# Constraints: satisfy demand for each piece size
for i, size in enumerate(sizes):
    prob += pulp.lpSum(pattern[i] * x[j] for j, pattern in enumerate(patterns)) >= demand[size]

# Solve
prob.solve()

# Collect and display the solution
used_patterns = [(j, x[j].varValue, patterns[j]) for j in range(len(patterns)) if x[j].varValue > 0]
used_patterns = sorted(used_patterns, key=lambda tup: -tup[1])

print("Optimal number of stock rods used:", pulp.value(prob.objective))
print("\nCutting patterns used:")
for j, count, pattern in used_patterns:
    print(f"Use pattern {pattern} -> {int(count)} times")

🧠 Code Explanation

generate_patterns: creates all combinations of pieces $(20, 45, 50)$ that can fit into a $100$-unit rod.
pulp.LpVariable: defines the variables $( x_j )$, how many times each pattern is used.
Constraints: ensure each type of piece meets its demand.
Objective: minimize the total number of rods used.

We now visualize the solution.

📊 Visualizing the Solution

Let’s create a simple bar chart to show how each rod is used and what pieces are cut from it.

# Visualization
colors = ['#ff9999','#66b3ff','#99ff99']
labels = [f"{s} unit" for s in sizes]

fig, ax = plt.subplots(figsize=(10, 6))

for i, (j, count, pattern) in enumerate(used_patterns):
    for k in range(int(count)):
        y = i + k * 0.1
        start = 0
        for size_idx, num in enumerate(pattern):
            for _ in range(num):
                ax.barh(y, sizes[size_idx], left=start, height=0.08, color=colors[size_idx])
                start += sizes[size_idx]

ax.set_title("Cutting Stock Solution Visualization")
ax.set_xlabel("Rod Length (100 units)")
ax.set_ylabel("Pattern Instances")
ax.set_xlim(0, stock_length)
ax.legend(labels, loc='upper right')
plt.tight_layout()
plt.show()

📈 Results

Optimal number of stock rods used: 40.0

Cutting patterns used:
Use pattern (0, 2, 0) -> 17 times
Use pattern (0, 0, 2) -> 12 times
Use pattern (5, 0, 0) -> 9 times
Use pattern (2, 1, 0) -> 2 times

📌 Interpretation of Results

The bar chart shows how each rod is used—what lengths are cut and how many rods follow the same pattern.
The objective value gives the minimal number of $100$-unit rods required.
We now know the optimal cutting strategy that minimizes waste and satisfies all customer demands.

✅ Conclusion

The Cutting Stock Problem is a fantastic application of linear programming.
Using Python and $PuLP$, we efficiently solved a real-world problem and even visualized the solution.

This model can be expanded further using column generation techniques when the number of patterns is too large.
But for many practical cases, this brute-force pattern generation combined with $PuLP$ is good enough.

Optimizing Inventory Management with EOQ in Python

April 4, 2025

Solving an Inventory Management Problem with Python

Inventory management is a crucial aspect of business operations.
Companies must balance between having enough stock to meet demand while minimizing holding costs.
Today, we’ll look at a specific inventory management problem, solve it using Python, and visualize the results for a better understanding.

Problem Statement

A retail store sells a particular product, and we want to determine the optimal inventory policy using the Economic Order Quantity (EOQ) model.
The EOQ formula helps find the order quantity that minimizes total inventory costs, including ordering and holding costs.

The formula for EOQ is given by:

$$
EOQ = \sqrt{\frac{2DS}{H}}
$$

Where:

$( D )$ = Annual demand (units)
$( S )$ = Ordering cost per order ($)
$( H )$ = Holding cost per unit per year ($)

Given:

$( D = 10,000 )$ units
$( S = 50 )$ dollars per order
$( H = 2 )$ dollars per unit per year

We will calculate the EOQ, determine the total cost, and visualize how the order quantity affects costs.

Python Implementation

Let’s write a Python script to compute EOQ and analyze costs.

import numpy as np
import matplotlib.pyplot as plt

def calculate_eoq(D, S, H):
    return np.sqrt((2 * D * S) / H)

def total_cost(D, S, H, Q):
    ordering_cost = (D / Q) * S
    holding_cost = (Q / 2) * H
    return ordering_cost + holding_cost

# Given data
D = 10000  # Annual demand
S = 50     # Ordering cost per order
H = 2      # Holding cost per unit per year

# Compute EOQ
eoq = calculate_eoq(D, S, H)
cost_at_eoq = total_cost(D, S, H, eoq)

print(f"Economic Order Quantity (EOQ): {eoq:.2f} units")
print(f"Total cost at EOQ: ${cost_at_eoq:.2f}")

# Visualizing the cost function
Q_values = np.linspace(50, 1000, 100)  # Order quantity range
costs = [total_cost(D, S, H, Q) for Q in Q_values]

plt.figure(figsize=(8, 5))
plt.plot(Q_values, costs, label='Total Cost', color='b')
plt.axvline(eoq, color='r', linestyle='--', label=f'EOQ: {eoq:.2f}')
plt.xlabel("Order Quantity (Q)")
plt.ylabel("Total Cost ($)")
plt.title("Total Inventory Cost vs Order Quantity")
plt.legend()
plt.grid()
plt.show()

Code Explanation

Calculate EOQ
- We use the EOQ formula to compute the optimal order quantity.
Calculate Total Cost
- The total cost function is derived from the sum of ordering and holding costs:
  - Ordering cost: $(\frac{D}{Q} S)$
  - Holding cost: $(\frac{Q}{2} H)$
- We compute this cost across different order quantities.
Visualization
- We plot how the total inventory cost varies with order quantity.
- The EOQ is marked with a red dashed line, showing the minimum cost point.

Results

Economic Order Quantity (EOQ): 707.11 units
Total cost at EOQ: $1414.21

Interpreting the Results

The EOQ value tells us how many units to order each time to minimize total cost.
The graph helps visualize why this is the optimal point—costs are lowest at EOQ.
If we order too frequently (small Q), ordering costs dominate.
If we order too infrequently (large Q), holding costs become excessive.

This analysis helps businesses make data-driven inventory decisions, optimizing costs and improving efficiency!

Optimizing Transportation Routes with Python and Linear Programming

April 3, 2025

Introduction

Transportation route optimization is a crucial problem in logistics and supply chain management.
The goal is to determine the most cost-effective way to transport goods from suppliers to customers while minimizing costs.
This type of problem is often formulated as a linear programming (LP) problem and solved using optimization techniques.

In this blog post, we will go through a practical example of solving a transportation problem using Python in Google Colaboratory.
We will utilize the PuLP library for linear programming and visualize the results using Matplotlib.

Problem Statement

A company needs to transport goods from two warehouses (suppliers) to three retail stores (destinations).
The transportation cost per unit of goods between each supplier and store is given in the table below.

Cost Table (per unit of goods)

	Store A	Store B	Store C
Warehouse X	$2	$3	$1
Warehouse Y	$5	$4	$8

Supply & Demand

Warehouse X can supply 50 units.
Warehouse Y can supply 60 units.
Store A requires 30 units.
Store B requires 40 units.
Store C requires 40 units.

The objective is to minimize the total transportation cost while ensuring that supply and demand constraints are met.

Mathematical Formulation

Define $ x_{ij} $ as the number of units transported from warehouse $ i $ to store $ j $.

Objective Function

We aim to minimize the total transportation cost:

$$
\text{Minimize } Z = 2x_{XA} + 3x_{XB} + 1x_{XC} + 5x_{YA} + 4x_{YB} + 8x_{YC}
$$

Constraints

Supply constraints (warehouses):
$$
x_{XA} + x_{XB} + x_{XC} \leq 50
$$
$$
x_{YA} + x_{YB} + x_{YC} \leq 60
$$

Demand constraints (stores):
$$
x_{XA} + x_{YA} = 30
$$
$$
x_{XB} + x_{YB} = 40
$$
$$
x_{XC} + x_{YC} = 40
$$

Non-negativity constraint:
$$
x_{ij} \geq 0
$$

Solving with Python (PuLP)

Let’s now solve this problem using Python.
Below is the complete code that sets up and solves the transportation optimization problem.

from pulp import LpProblem, LpMinimize, LpVariable, lpSum
import matplotlib.pyplot as plt
import numpy as np

# Define the problem
problem = LpProblem("Transportation_Optimization", LpMinimize)

# Define decision variables
x_XA = LpVariable("x_XA", lowBound=0, cat='Continuous')
x_XB = LpVariable("x_XB", lowBound=0, cat='Continuous')
x_XC = LpVariable("x_XC", lowBound=0, cat='Continuous')
x_YA = LpVariable("x_YA", lowBound=0, cat='Continuous')
x_YB = LpVariable("x_YB", lowBound=0, cat='Continuous')
x_YC = LpVariable("x_YC", lowBound=0, cat='Continuous')

# Objective function: Minimize total cost
problem += (2*x_XA + 3*x_XB + 1*x_XC + 
            5*x_YA + 4*x_YB + 8*x_YC), "Total Transportation Cost"

# Supply constraints
problem += x_XA + x_XB + x_XC <= 50, "Supply_X"
problem += x_YA + x_YB + x_YC <= 60, "Supply_Y"

# Demand constraints
problem += x_XA + x_YA == 30, "Demand_A"
problem += x_XB + x_YB == 40, "Demand_B"
problem += x_XC + x_YC == 40, "Demand_C"

# Solve the problem
problem.solve()

# Output results
print("\nOptimal Transportation Plan:")
for v in problem.variables():
    print(f"{v.name} = {v.varValue}")
print(f"\nTotal Minimum Cost: ${problem.objective.value()}")

# Visualization
labels = ["X → A", "X → B", "X → C", "Y → A", "Y → B", "Y → C"]
values = [v.varValue for v in problem.variables()]

plt.figure(figsize=(8, 6))
plt.bar(labels, values, color=['blue', 'green', 'red', 'purple', 'orange', 'brown'])
plt.xlabel("Routes")
plt.ylabel("Units Transported")
plt.title("Optimized Transportation Plan")
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.show()

Explanation of the Code

Defining the Problem
- We create a linear programming problem using LpProblem().
- The goal is minimization of transportation cost.
Defining Decision Variables
- Each route (e.g., x_XA, x_XB, etc.) is a continuous variable representing the number of units transported.
Setting the Objective Function
- We minimize the total cost by summing cost × transported units for each route.
Defining Constraints
- Supply constraints ensure warehouses do not exceed their limits.
- Demand constraints ensure stores receive the required amount.
Solving the Problem
- The solve() method finds the optimal transportation plan.
Output Results & Visualization
- The results are printed in a readable format.
- A bar chart is generated using Matplotlib to visualize the transportation flow.

Results & Interpretation

After running the script, we obtain:

Optimal Transportation Plan:
x_XA = 10.0
x_XB = 0.0
x_XC = 40.0
x_YA = 20.0
x_YB = 40.0
x_YC = 0.0

Total Minimum Cost: $320.0

Key Insights

Warehouse X supplies 10 units to Store A and 40 units to Store C.
Warehouse Y supplies 20 units to Store A and 40 units to Store B.
Total cost is minimized to $320.0.

Graph Interpretation

The bar chart shows how many units each route transports.
The highest transportation occurs between X → C and Y → B.
The optimized plan avoids unnecessary high-cost routes.

Conclusion

This example demonstrates how linear programming effectively optimizes transportation routes in logistics.
By using Python’s PuLP, we minimized costs while satisfying all constraints.

Optimizing Operating Room Scheduling with Python: A Practical Approach

🏥 The Problem: Operating Room Scheduling

📦 Data Setup and Model in Python

✅ Explanation of the Code

📊 Visualizing the Results

📈 Interpretation of the Output

🧠 Takeaways

🆘 Optimizing Evacuation Routes and Resource Allocation During Disasters with Python

🧠 The Problem: Evacuation and Supply Allocation

🧮 Mathematical Formulation

Objective:

Subject to:

🧪 Let’s Code!

✅ Install PuLP

🧮 Python Code (Evacuation Optimization)

📊 Visualization

🔍 Results and Analysis

🧯 Next Steps

🛠️ Designing a Reliable Power Grid with Python

📘 The Problem Statement

🧮 Mathematical Formulation

🐍 Python Code

🔍 Code Breakdown

📊 Results & Interpretation

✅ Conclusion

Solving the Capital Allocation Problem with Python

Problem Statement

Mathematical Formulation

Code Explanation

Results and Visualization

Insights and Conclusions

🛠️ Solving the Machine Assignment Problem with Python: A Step-by-Step Guide

📘 Problem Statement

🔍 Mathematical Formulation

🧠 Python Implementation

🧩 Code Explanation

📊 Visualizing the Assignment

📌 Result Interpretation

✅ Conclusion

Understanding the Vehicle Routing Problem

A Practical VRP Example

Code Explanation

1. VRPSolver Class

2. Solution Approaches

3. Visualization

Results Analysis

Optimization Opportunities

Conclusion

🧠 Solving a Network Design Problem with Python

📌 Problem Statement: Minimum Cost Network Design

📊 Sample Data

🔢 Mathematical Formulation

🐍 Python Code: Prim’s Algorithm

🧠 Code Explanation

📈 Visualizing the Network

📌 Final Result

🧩 Conclusion

Solving the Cutting Stock Problem with Python 🧠✂️

🧩 Problem Statement

📌 Mathematical Formulation

🐍 Python Implementation

🧠 Code Explanation

📊 Visualizing the Solution

📈 Results

📌 Interpretation of Results

✅ Conclusion

Solving an Inventory Management Problem with Python

Problem Statement

Python Implementation

Code Explanation

Results

Interpreting the Results

Introduction

Problem Statement

Cost Table (per unit of goods)

Supply & Demand

Mathematical Formulation

Objective Function

Constraints

Solving with Python (PuLP)

1. `VRPSolver` Class