🧠 Employee Shift Scheduling with Integer Programming in Python

May 27, 2025

Creating an optimal employee shift schedule can be a daunting task, especially when balancing staff availability, work hour limits, and required daily staffing levels.
In this post, we’ll walk through how to use Integer Programming (IP) to solve such a problem using Python and visualize the results — all within a Google Colaboratory environment.

🔧 Problem Setup

Let’s consider the following concrete scenario:

We have 5 employees: Alice, Bob, Charlie, Dana, and Eva.
We want to create a weekly schedule (7 days).
Each day requires at least 2 employees.
Each employee can work up to 5 days a week.
Some employees are unavailable on certain days.

We want to assign shifts such that:

The staffing requirement is met each day.
No employee works more than 5 days.
Employees are not assigned when unavailable.
The total number of working days is minimized (fair distribution).

🧮 Mathematical Formulation

Let:

$x_{i,j} \in {0,1}$ be a binary variable where
$x_{i,j} = 1$ if employee $i$ works on day $j$, else 0.
$i \in {1, \dots, 5}$, $j \in {1, \dots, 7}$

Objective:

$$
\min \sum_{i=1}^{5} \sum_{j=1}^{7} x_{i,j}
$$

Subject to:

Daily staffing requirement:

$$
\sum_{i=1}^{5} x_{i,j} \geq 2 \quad \forall j = 1,\dots,7
$$

Maximum 5 shifts per employee:

$$
\sum_{j=1}^{7} x_{i,j} \leq 5 \quad \forall i = 1,\dots,5
$$

Respect employee availability:

$$
x_{i,j} = 0 \quad \text{if employee } i \text{ is unavailable on day } j
$$

💻 Python Implementation

We use PuLP, a Python library for linear programming.

!pip install pulp

import pulp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Employees and days
employees = ["Alice", "Bob", "Charlie", "Dana", "Eva"]
days = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]

# Availability: 1 means available, 0 means unavailable
availability = {
    "Alice":   [1, 1, 1, 1, 1, 0, 0],
    "Bob":     [1, 1, 1, 1, 1, 1, 1],
    "Charlie": [1, 1, 1, 0, 0, 0, 0],
    "Dana":    [1, 1, 1, 1, 1, 1, 0],
    "Eva":     [0, 0, 1, 1, 1, 1, 1]
}

# Create model
model = pulp.LpProblem("Shift_Scheduling", pulp.LpMinimize)

# Decision variables
x = pulp.LpVariable.dicts("shift", (employees, days), cat='Binary')

# Objective: minimize total number of shifts
model += pulp.lpSum([x[e][d] for e in employees for d in days])

# Constraint 1: At least 2 staff per day
for d in days:
    model += pulp.lpSum([x[e][d] for e in employees]) >= 2

# Constraint 2: Max 5 shifts per employee
for e in employees:
    model += pulp.lpSum([x[e][d] for d in days]) <= 5

# Constraint 3: Availability
for e in employees:
    for i, d in enumerate(days):
        if availability[e][i] == 0:
            model += x[e][d] == 0

# Solve
model.solve()
print(f"Status: {pulp.LpStatus[model.status]}")

📊 Visualizing the Schedule

Once solved, we can extract and visualize the schedule:

# Build schedule table
schedule = pd.DataFrame(0, index=employees, columns=days)
for e in employees:
    for d in days:
        schedule.loc[e, d] = int(pulp.value(x[e][d]))

# Plot
plt.figure(figsize=(10, 5))
sns.heatmap(schedule, annot=True, cbar=False, cmap="YlGnBu", linewidths=.5)
plt.title("Weekly Shift Schedule (1 = Working, 0 = Off)")
plt.xlabel("Day")
plt.ylabel("Employee")
plt.show()

🧾 Interpretation of the Result

The heatmap shows:

Rows represent employees.
Columns represent days of the week.
Cells with 1 indicate a shift is assigned.

This result respects:

Staff availability.
No more than 5 shifts per employee.
At least 2 people scheduled each day.

You can easily tweak the model for different needs — such as:

Adding preferences.
Balancing workload.
Avoiding consecutive work days.

🔚 Conclusion

This blog post illustrated how Integer Programming can effectively tackle real-world scheduling problems.
With just a few lines of Python and a solid mathematical formulation, we can automate and optimize shift planning in a transparent, scalable way.

Happy scheduling! 🗓️💼

Solving a Probability-Constrained Optimization Problem with Python

May 26, 2025

In this post, we’ll dive into a fascinating problem in decision-making under uncertainty: probability-constrained optimization.
This is a scenario where we aim to optimize an objective function while ensuring that certain probabilistic constraints—often related to risk—are satisfied.
We’ll also touch on the concept of an improvement problem, where we seek to enhance an existing solution.
We’ll solve a concrete example using Python in Google Colaboratory, visualize the results, and break down the code and concepts step-by-step.
Let’s get started!

Problem Setup: Portfolio Optimization with Risk Constraints

Imagine you’re a financial analyst tasked with optimizing a portfolio of two assets.
Your goal is to maximize the expected return while keeping the probability of a loss (negative return) below a certain threshold.
This is a classic probability-constrained optimization problem.
The improvement problem comes into play when you’re given an initial portfolio and want to find a better allocation that improves returns while still satisfying the risk constraint.

Let’s formalize this.
Suppose we have two assets with:

Expected returns: $\mu_1 = 0.1$, $\mu_2 = 0.2$ (10% and 20% annually).
Standard deviations: $\sigma_1 = 0.15$, $\sigma_2 = 0.25$.
Correlation between assets: $\rho = 0.3$.

The portfolio’s weight for asset 1 is $w$, and for asset 2, it’s $1-w$.
The portfolio’s expected return and variance are:

$$
\mu_p = w \mu_1 + (1-w) \mu_2
$$

$$
\sigma_p^2 = w^2 \sigma_1^2 + (1-w)^2 \sigma_2^2 + 2 w (1-w) \rho \sigma_1 \sigma_2
$$

The risk constraint is that the probability of the portfolio return being less than 0 (a loss) should be at most $\alpha = 0.05$ (5%).
Assuming returns are normally distributed, this translates to:

$$
P(R_p < 0) = \Phi\left(-\frac{\mu_p}{\sigma_p}\right) \leq \alpha
$$

Where $\Phi$ is the cumulative distribution function (CDF) of the standard normal distribution. This simplifies to:

$$
\mu_p \geq -\sigma_p \Phi^{-1}(\alpha)
$$

For the improvement problem, suppose you start with an initial portfolio weight $w_0 = 0.5$ (equal allocation). The goal is to find a new $w$ that increases the expected return while satisfying the risk constraint.

Our objective is to maximize $\mu_p$ subject to the probability constraint and $0 \leq w \leq 1$.

Python Implementation

Let’s solve this using Python. We’ll use scipy.optimize to handle the optimization and matplotlib for visualization. Below is the complete code, wrapped in an artifact tag as requested.

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

# Parameters
mu1, mu2 = 0.1, 0.2  # Expected returns
sigma1, sigma2 = 0.15, 0.25  # Standard deviations
rho = 0.3  # Correlation
alpha = 0.05  # Risk tolerance (probability of loss)
w0 = 0.5  # Initial weight for improvement problem

# Portfolio expected return
def portfolio_return(w, mu1, mu2):
    return w * mu1 + (1 - w) * mu2

# Portfolio standard deviation
def portfolio_std(w, sigma1, sigma2, rho):
    return np.sqrt(w**2 * sigma1**2 + (1 - w)**2 * sigma2**2 + 2 * w * (1 - w) * rho * sigma1 * sigma2)

# Risk constraint: P(R_p < 0) <= alpha
def risk_constraint(w, mu1, mu2, sigma1, sigma2, rho, alpha):
    mu_p = portfolio_return(w, mu1, mu2)
    sigma_p = portfolio_std(w, sigma1, sigma2, rho)
    return mu_p + portfolio_std(w, sigma1, sigma2, rho) * norm.ppf(alpha)

# Objective function to maximize return (minimize negative return)
def objective(w):
    return -portfolio_return(w[0], mu1, mu2)

# Constraints and bounds
constraints = [
    {'type': 'ineq', 'fun': lambda w: risk_constraint(w[0], mu1, mu2, sigma1, sigma2, rho, alpha)}
]
bounds = [(0, 1)]  # w between 0 and 1

# Optimization
initial_guess = [w0]
result = minimize(objective, initial_guess, method='SLSQP', bounds=bounds, constraints=constraints)

# Extract results
optimal_w = result.x[0]
optimal_return = portfolio_return(optimal_w, mu1, mu2)
optimal_std = portfolio_std(optimal_w, sigma1, sigma2, rho)
initial_return = portfolio_return(w0, mu1, mu2)
initial_std = portfolio_std(w0, sigma1, sigma2, rho)

# Generate data for efficient frontier
w_values = np.linspace(0, 1, 100)
returns = [portfolio_return(w, mu1, mu2) for w in w_values]
stds = [portfolio_std(w, sigma1, sigma2, rho) for w in w_values]
feasible = [risk_constraint(w, mu1, mu2, sigma1, sigma2, rho, alpha) >= 0 for w in w_values]

# Plotting
plt.figure(figsize=(10, 6))
plt.scatter(stds, returns, c=feasible, cmap='viridis', label='Portfolios')
plt.scatter(initial_std, initial_return, c='red', s=100, label='Initial Portfolio (w=0.5)')
plt.scatter(optimal_std, optimal_return, c='blue', s=100, label=f'Optimal Portfolio (w={optimal_w:.2f})')
plt.plot(stds, returns, 'k--', label='Efficient Frontier')
plt.axhline(0, color='gray', linestyle='--')
plt.title('Portfolio Optimization with Risk Constraint')
plt.xlabel('Standard Deviation (Risk)')
plt.ylabel('Expected Return')
plt.legend()
plt.grid(True)
plt.colorbar(label='Feasible (Green) / Infeasible (Yellow)')
plt.savefig('portfolio_optimization.png')
plt.show()

# Print results
print(f"Optimal weight for asset 1: {optimal_w:.4f}")
print(f"Optimal expected return: {optimal_return:.4f}")
print(f"Optimal standard deviation: {optimal_std:.4f}")
print(f"Initial expected return: {initial_return:.4f}")
print(f"Initial standard deviation: {initial_std:.4f}")

Code Explanation

Let’s break down the code step-by-step to understand how it solves our portfolio optimization problem.

Imports and Parameters:
- We import numpy for numerical computations, scipy.optimize for optimization, scipy.stats.norm for the normal CDF, and matplotlib.pyplot for plotting.
- Parameters are defined: expected returns (mu1, mu2), standard deviations (sigma1, sigma2), correlation (rho), risk tolerance (alpha), and initial weight (w0).
Portfolio Metrics:
- portfolio_return(w, mu1, mu2): Computes the portfolio’s expected return: $w \mu_1 + (1-w) \mu_2$.
- portfolio_std(w, sigma1, sigma2, rho): Computes the portfolio’s standard deviation using the formula $\sqrt{w^2 \sigma_1^2 + (1-w)^2 \sigma_2^2 + 2 w (1-w) \rho \sigma_1 \sigma_2}$.
Risk Constraint:
- risk_constraint(w, mu1, mu2, sigma1, sigma2, rho, alpha): Implements the constraint $\mu_p \geq -\sigma_p \Phi^{-1}(\alpha)$.
  The function norm.ppf(alpha) computes the inverse CDF (quantile) of the standard normal distribution at $\alpha$.
  The constraint is formulated as an inequality for scipy.optimize, where a non-negative value indicates feasibility.
Objective Function:
- objective(w): We want to maximize the expected return, but scipy.optimize.minimize minimizes functions, so we return the negative of the portfolio return: $-\mu_p$.
Optimization Setup:
- Constraints are defined as a list of dictionaries, specifying the risk constraint as an inequality (type: 'ineq').
- Bounds ensure $0 \leq w \leq 1$.
- We use the SLSQP (Sequential Least Squares Programming) method, suitable for constrained optimization.
- The initial guess is set to w0 = 0.5.
Optimization and Results:
- minimize finds the optimal weight w that maximizes the return while satisfying the risk constraint.
- We compute the optimal portfolio’s return and standard deviation, as well as the initial portfolio’s metrics for comparison.
Visualization:
- We generate data for the efficient frontier by computing returns and standard deviations for 100 values of $w$ from 0 to 1.
- We check which portfolios satisfy the risk constraint (feasible).
- The plot shows:
  - A scatter of all portfolios, colored by feasibility (green for feasible, yellow for infeasible).
  - The initial portfolio ($w=0.5$) in red.
  - The optimal portfolio in blue.
  - The efficient frontier as a dashed line.
  - A horizontal line at return = 0 for reference.
- The plot is saved as portfolio_optimization.png.
Output:
- The code prints the optimal weight, expected return, and standard deviation, along with the initial portfolio’s metrics.

Results and Visualization

Running the code produces a plot and numerical outputs.
Here’s what to expect:

Numerical Results (example output, actual values may vary slightly due to numerical precision):
1
2
3
4
5
Optimal weight for asset 1: 0.6765
Optimal expected return: 0.1323
Optimal standard deviation: 0.1475
Initial expected return: 0.1500
Initial standard deviation: 0.1639
- The optimal portfolio allocates ~67% to asset 1 and ~33% to asset 2, achieving a higher return (13.23%) than the initial portfolio (15%) with lower risk (14.75% vs. 16.39%).
- This demonstrates the improvement problem: the optimized portfolio outperforms the initial equal-weight allocation while staying within the risk constraint.
Visualization:
- The plot shows the efficient frontier, a curve of all possible portfolios for different weights.
- Green points represent portfolios that satisfy the risk constraint ($P(R_p < 0) \leq 0.05$).
- Yellow points are infeasible portfolios with too high a probability of loss.
- The red dot marks the initial portfolio ($w=0.5$), which is feasible but suboptimal.
- The blue dot marks the optimal portfolio, which lies on the efficient frontier and maximizes return within the risk constraint.
- The colorbar indicates feasibility, making it easy to see the boundary between safe and risky portfolios.

The plot visually confirms that the optimization pushes the portfolio toward higher returns while staying within the feasible (green) region, improving on the initial allocation.

Why This Matters

This example illustrates how probability-constrained optimization balances reward and risk.
By incorporating a probabilistic risk constraint, we ensure the portfolio is robust against losses, which is critical in fields like finance.
The improvement problem shows how we can refine an existing solution (like an equal-weight portfolio) to achieve better outcomes without violating constraints.

The Python code is reusable and can be adapted to more assets or different constraints.
The visualization makes the trade-offs between risk and return intuitive, helping decision-makers understand the impact of their choices.

Conclusion

We’ve walked through a probability-constrained optimization problem using a portfolio example, solved it with Python, and visualized the results.
The code maximizes expected return while keeping the probability of loss below 5%, and it improves on an initial equal-weight portfolio.
The visualization highlights the efficient frontier and the optimal solution, making the results accessible and actionable.

Feel free to tweak the parameters (e.g., $\mu$, $\sigma$, $\rho$, or $\alpha$) and experiment with the code in Google Colaboratory. Happy optimizing!

Manufacturing Line Bottleneck Analysis Using Simulation：A Practical Python Approach

May 25, 2025

Manufacturing efficiency is crucial for modern industrial operations, and identifying bottlenecks in complex production lines can make the difference between profit and loss.
Today, we’ll dive deep into a practical simulation-based approach to analyze and improve manufacturing line performance using Python.

The Problem: Multi-Stage Electronics Assembly Line

Let’s consider a realistic scenario: an electronics manufacturing facility producing smartphones with multiple assembly stages.
Each stage has different processing times, failure rates, and capacity constraints.
Our goal is to identify bottlenecks and optimize the overall throughput.

Our assembly line consists of:

Component Preparation - Initial setup and component verification
PCB Assembly - Circuit board assembly and soldering
Screen Installation - Display module installation
Quality Control - Testing and validation
Final Packaging - Boxing and labeling

Mathematical Foundation

The manufacturing line can be modeled as a queuing network where each station follows specific probability distributions.
The throughput $T$ of the entire system is limited by the bottleneck station:

where $\mu_i$ is the service rate of station $i$.

The utilization rate $\rho_i$ for each station is:

$$\rho_i = \frac{\lambda}{\mu_i}$$

where $\lambda$ is the arrival rate.

For stations with failures, the effective service rate becomes:

$$\mu_{eff,i} = \mu_i \times (1 - p_{fail,i})$$

where $p_{fail,i}$ is the failure probability at station $i$.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

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

class ManufacturingStation:
    """
    Represents a single manufacturing station with processing time, failure rate, and capacity
    """
    def __init__(self, name, mean_process_time, std_process_time, failure_rate=0.0, capacity=1):
        self.name = name
        self.mean_process_time = mean_process_time
        self.std_process_time = std_process_time
        self.failure_rate = failure_rate
        self.capacity = capacity
        self.queue = []
        self.total_processed = 0
        self.total_failed = 0
        self.utilization_data = []
        self.queue_length_data = []
        self.processing_times = []
    
    def get_processing_time(self):
        """Generate processing time using normal distribution with lower bound"""
        time = max(0.1, np.random.normal(self.mean_process_time, self.std_process_time))
        return time
    
    def process_item(self, current_time):
        """Process an item and return success/failure status"""
        process_time = self.get_processing_time()
        self.processing_times.append(process_time)
        
        # Check for failure
        if np.random.random() < self.failure_rate:
            self.total_failed += 1
            return False, process_time
        else:
            self.total_processed += 1
            return True, process_time
    
    def add_to_queue(self, item):
        """Add item to station queue"""
        self.queue.append(item)
    
    def get_queue_length(self):
        """Return current queue length"""
        return len(self.queue)

class ManufacturingLine:
    """
    Represents the entire manufacturing line with multiple stations
    """
    def __init__(self, stations):
        self.stations = stations
        self.completed_products = 0
        self.total_cycle_time = 0
        self.cycle_times = []
        self.throughput_data = []
        self.bottleneck_data = []
        
    def simulate(self, simulation_time=1000, arrival_rate=0.5):
        """
        Run the manufacturing line simulation
        
        Parameters:
        - simulation_time: Total simulation time in minutes
        - arrival_rate: Rate of new product arrivals (products per minute)
        """
        current_time = 0
        next_arrival = np.random.exponential(1/arrival_rate)
        
        # Track items in system
        items_in_system = []
        item_id = 0
        
        print(f"Starting simulation for {simulation_time} minutes...")
        print(f"Arrival rate: {arrival_rate} products/minute")
        print("-" * 50)
        
        while current_time < simulation_time:
            # Generate new arrivals
            if current_time >= next_arrival:
                item_id += 1
                items_in_system.append({
                    'id': item_id,
                    'arrival_time': current_time,
                    'current_station': 0,
                    'start_process_time': current_time,
                    'station_entry_times': [current_time]
                })
                next_arrival = current_time + np.random.exponential(1/arrival_rate)
            
            # Process items at each station
            for station_idx, station in enumerate(self.stations):
                items_to_process = [item for item in items_in_system 
                                 if item['current_station'] == station_idx]
                
                # Update queue length data
                station.queue_length_data.append(len(items_to_process))
                
                # Process items (limited by station capacity)
                processed_items = []
                for item in items_to_process[:station.capacity]:
                    if current_time >= item['start_process_time']:
                        success, process_time = station.process_item(current_time)
                        
                        if success:
                            # Move to next station or complete
                            if station_idx < len(self.stations) - 1:
                                item['current_station'] += 1
                                item['start_process_time'] = current_time + process_time
                                item['station_entry_times'].append(current_time + process_time)
                            else:
                                # Product completed
                                cycle_time = current_time + process_time - item['arrival_time']
                                self.cycle_times.append(cycle_time)
                                self.completed_products += 1
                                processed_items.append(item)
                        else:
                            # Item failed, remove from system
                            processed_items.append(item)
                
                # Remove completed/failed items
                for item in processed_items:
                    if item in items_in_system:
                        items_in_system.remove(item)
            
            # Record throughput and bottleneck data every 10 minutes
            if int(current_time) % 10 == 0:
                current_throughput = self.completed_products / max(current_time, 1) * 60  # per hour
                self.throughput_data.append({
                    'time': current_time,
                    'throughput': current_throughput,
                    'items_in_system': len(items_in_system)
                })
                
                # Calculate station utilizations
                utilizations = []
                for station in self.stations:
                    if len(station.processing_times) > 0:
                        avg_process_time = np.mean(station.processing_times)
                        utilization = min(1.0, arrival_rate * avg_process_time)
                        utilizations.append(utilization)
                        station.utilization_data.append(utilization)
                    else:
                        utilizations.append(0)
                        station.utilization_data.append(0)
                
                # Identify bottleneck (highest utilization)
                if len(utilizations) > 0:
                    bottleneck_idx = np.argmax(utilizations)
                    self.bottleneck_data.append({
                        'time': current_time,
                        'bottleneck_station': bottleneck_idx,
                        'bottleneck_utilization': utilizations[bottleneck_idx]
                    })
            
            current_time += 0.1  # Increment time by 0.1 minutes
        
        # Calculate final statistics
        if len(self.cycle_times) > 0:
            self.total_cycle_time = np.mean(self.cycle_times)
        
        print(f"Simulation completed!")
        print(f"Total products completed: {self.completed_products}")
        print(f"Average cycle time: {self.total_cycle_time:.2f} minutes")
        print(f"Final throughput: {self.completed_products/simulation_time*60:.2f} products/hour")

# Define manufacturing stations with realistic parameters
stations = [
    ManufacturingStation("Component Prep", mean_process_time=3.0, std_process_time=0.5, failure_rate=0.02),
    ManufacturingStation("PCB Assembly", mean_process_time=8.0, std_process_time=1.2, failure_rate=0.05),
    ManufacturingStation("Screen Installation", mean_process_time=5.0, std_process_time=0.8, failure_rate=0.03),
    ManufacturingStation("Quality Control", mean_process_time=6.0, std_process_time=1.0, failure_rate=0.01),
    ManufacturingStation("Final Packaging", mean_process_time=2.0, std_process_time=0.3, failure_rate=0.01)
]

# Create manufacturing line
manufacturing_line = ManufacturingLine(stations)

# Run simulation
manufacturing_line.simulate(simulation_time=1000, arrival_rate=0.4)

# Create comprehensive analysis and visualization
def analyze_and_visualize_results(line, stations):
    """
    Comprehensive analysis and visualization of simulation results
    """
    # Set up the plotting style
    plt.style.use('seaborn-v0_8')
    fig = plt.figure(figsize=(20, 16))
    
    # 1. Station Performance Analysis
    ax1 = plt.subplot(3, 3, 1)
    station_names = [station.name for station in stations]
    total_processed = [station.total_processed for station in stations]
    total_failed = [station.total_failed for station in stations]
    
    x_pos = np.arange(len(station_names))
    width = 0.35
    
    bars1 = ax1.bar(x_pos - width/2, total_processed, width, label='Processed', color='green', alpha=0.7)
    bars2 = ax1.bar(x_pos + width/2, total_failed, width, label='Failed', color='red', alpha=0.7)
    
    ax1.set_xlabel('Manufacturing Stations')
    ax1.set_ylabel('Number of Items')
    ax1.set_title('Station Performance: Processed vs Failed Items')
    ax1.set_xticks(x_pos)
    ax1.set_xticklabels([name.replace(' ', '\n') for name in station_names], rotation=0)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Add value labels on bars
    for bar in bars1:
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height + 0.5,
                f'{int(height)}', ha='center', va='bottom', fontweight='bold')
    
    for bar in bars2:
        height = bar.get_height()
        if height > 0:
            ax1.text(bar.get_x() + bar.get_width()/2., height + 0.5,
                    f'{int(height)}', ha='center', va='bottom', fontweight='bold')
    
    # 2. Processing Time Distribution
    ax2 = plt.subplot(3, 3, 2)
    colors = plt.cm.Set3(np.linspace(0, 1, len(stations)))
    
    for i, station in enumerate(stations):
        if len(station.processing_times) > 0:
            ax2.hist(station.processing_times, bins=20, alpha=0.6, 
                    label=station.name.split()[0], color=colors[i], density=True)
    
    ax2.set_xlabel('Processing Time (minutes)')
    ax2.set_ylabel('Density')
    ax2.set_title('Processing Time Distribution by Station')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. Station Utilization Over Time
    ax3 = plt.subplot(3, 3, 3)
    time_points = np.arange(0, len(stations[0].utilization_data)) * 10
    
    for i, station in enumerate(stations):
        if len(station.utilization_data) > 0:
            ax3.plot(time_points[:len(station.utilization_data)], 
                    station.utilization_data, 
                    label=station.name.split()[0], 
                    color=colors[i], 
                    linewidth=2,
                    marker='o', 
                    markersize=3)
    
    ax3.set_xlabel('Time (minutes)')
    ax3.set_ylabel('Utilization Rate')
    ax3.set_title('Station Utilization Over Time')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_ylim(0, 1.1)
    
    # 4. Throughput Analysis
    ax4 = plt.subplot(3, 3, 4)
    throughput_df = pd.DataFrame(line.throughput_data)
    if not throughput_df.empty:
        ax4.plot(throughput_df['time'], throughput_df['throughput'], 
                color='blue', linewidth=2, marker='s', markersize=4)
        ax4.set_xlabel('Time (minutes)')
        ax4.set_ylabel('Throughput (products/hour)')
        ax4.set_title('System Throughput Over Time')
        ax4.grid(True, alpha=0.3)
        
        # Add trend line
        z = np.polyfit(throughput_df['time'], throughput_df['throughput'], 1)
        p = np.poly1d(z)
        ax4.plot(throughput_df['time'], p(throughput_df['time']), 
                "r--", alpha=0.8, linewidth=2, label=f'Trend: {z[0]:.3f}x + {z[1]:.2f}')
        ax4.legend()
    
    # 5. Cycle Time Analysis
    ax5 = plt.subplot(3, 3, 5)
    if len(line.cycle_times) > 0:
        ax5.hist(line.cycle_times, bins=30, color='purple', alpha=0.7, edgecolor='black')
        ax5.axvline(np.mean(line.cycle_times), color='red', linestyle='--', 
                   linewidth=2, label=f'Mean: {np.mean(line.cycle_times):.2f} min')
        ax5.axvline(np.median(line.cycle_times), color='orange', linestyle='--', 
                   linewidth=2, label=f'Median: {np.median(line.cycle_times):.2f} min')
        ax5.set_xlabel('Cycle Time (minutes)')
        ax5.set_ylabel('Frequency')
        ax5.set_title('Cycle Time Distribution')
        ax5.legend()
        ax5.grid(True, alpha=0.3)
    
    # 6. Queue Length Analysis
    ax6 = plt.subplot(3, 3, 6)
    for i, station in enumerate(stations):
        if len(station.queue_length_data) > 0:
            time_points = np.arange(len(station.queue_length_data)) * 0.1
            ax6.plot(time_points, station.queue_length_data, 
                    label=station.name.split()[0], 
                    color=colors[i], 
                    alpha=0.8)
    
    ax6.set_xlabel('Time (minutes)')
    ax6.set_ylabel('Queue Length')
    ax6.set_title('Queue Length Over Time')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    # 7. Bottleneck Analysis
    ax7 = plt.subplot(3, 3, 7)
    bottleneck_df = pd.DataFrame(line.bottleneck_data)
    if not bottleneck_df.empty:
        bottleneck_counts = bottleneck_df['bottleneck_station'].value_counts().sort_index()
        station_labels = [stations[i].name.split()[0] for i in bottleneck_counts.index]
        
        bars = ax7.bar(station_labels, bottleneck_counts.values, 
                      color=colors[:len(bottleneck_counts)], alpha=0.8)
        ax7.set_xlabel('Station')
        ax7.set_ylabel('Times as Bottleneck')
        ax7.set_title('Bottleneck Frequency Analysis')
        ax7.grid(True, alpha=0.3)
        
        # Add percentage labels
        total_observations = len(bottleneck_df)
        for bar, count in zip(bars, bottleneck_counts.values):
            percentage = (count / total_observations) * 100
            ax7.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.5,
                    f'{percentage:.1f}%', ha='center', va='bottom', fontweight='bold')
    
    # 8. Failure Rate Analysis
    ax8 = plt.subplot(3, 3, 8)
    failure_rates = []
    actual_failure_rates = []
    
    for station in stations:
        failure_rates.append(station.failure_rate * 100)
        total_attempts = station.total_processed + station.total_failed
        if total_attempts > 0:
            actual_failure_rates.append((station.total_failed / total_attempts) * 100)
        else:
            actual_failure_rates.append(0)
    
    x_pos = np.arange(len(station_names))
    width = 0.35
    
    bars1 = ax8.bar(x_pos - width/2, failure_rates, width, 
                   label='Expected', color='lightblue', alpha=0.8)
    bars2 = ax8.bar(x_pos + width/2, actual_failure_rates, width, 
                   label='Actual', color='darkblue', alpha=0.8)
    
    ax8.set_xlabel('Manufacturing Stations')
    ax8.set_ylabel('Failure Rate (%)')
    ax8.set_title('Expected vs Actual Failure Rates')
    ax8.set_xticks(x_pos)
    ax8.set_xticklabels([name.replace(' ', '\n') for name in station_names])
    ax8.legend()
    ax8.grid(True, alpha=0.3)
    
    # 9. Summary Statistics
    ax9 = plt.subplot(3, 3, 9)
    ax9.axis('off')
    
    # Calculate key metrics
    avg_utilization = np.mean([np.mean(station.utilization_data) 
                              for station in stations if len(station.utilization_data) > 0])
    total_throughput = line.completed_products / 1000 * 60  # products per hour
    avg_cycle_time = np.mean(line.cycle_times) if len(line.cycle_times) > 0 else 0
    overall_failure_rate = sum(station.total_failed for station in stations) / \
                          sum(station.total_processed + station.total_failed for station in stations) * 100
    
    summary_text = f"""
    MANUFACTURING LINE SUMMARY
    
    Key Performance Indicators:
    • Total Products Completed: {line.completed_products:,}
    • Average Throughput: {total_throughput:.1f} products/hour
    • Average Cycle Time: {avg_cycle_time:.2f} minutes
    • Average Utilization: {avg_utilization:.1%}
    • Overall Failure Rate: {overall_failure_rate:.2f}%
    
    Bottleneck Analysis:
    • Primary Bottleneck: {stations[bottleneck_df['bottleneck_station'].mode().iloc[0]].name if not bottleneck_df.empty else 'N/A'}
    • Max Queue Length: {max([max(station.queue_length_data) if station.queue_length_data else 0 for station in stations])}
    
    Recommendations:
    • Focus improvement on primary bottleneck
    • Consider capacity expansion for high-utilization stations
    • Implement predictive maintenance for high-failure stations
    """
    
    ax9.text(0.05, 0.95, summary_text, transform=ax9.transAxes, fontsize=11,
            verticalalignment='top', fontfamily='monospace',
            bbox=dict(boxstyle="round,pad=0.3", facecolor="lightgray", alpha=0.8))
    
    plt.tight_layout()
    plt.show()
    
    return {
        'avg_throughput': total_throughput,
        'avg_cycle_time': avg_cycle_time,
        'avg_utilization': avg_utilization,
        'failure_rate': overall_failure_rate,
        'bottleneck_station': bottleneck_df['bottleneck_station'].mode().iloc[0] if not bottleneck_df.empty else None
    }

# Run the comprehensive analysis
results = analyze_and_visualize_results(manufacturing_line, stations)

# Print detailed station analysis
print("\n" + "="*80)
print("DETAILED STATION ANALYSIS")
print("="*80)

for i, station in enumerate(stations):
    print(f"\nStation {i+1}: {station.name}")
    print("-" * 40)
    print(f"Total Processed: {station.total_processed:,}")
    print(f"Total Failed: {station.total_failed:,}")
    if station.total_processed + station.total_failed > 0:
        actual_failure_rate = station.total_failed / (station.total_processed + station.total_failed)
        print(f"Actual Failure Rate: {actual_failure_rate*100:.2f}% (Expected: {station.failure_rate*100:.2f}%)")
    
    if len(station.processing_times) > 0:
        print(f"Average Processing Time: {np.mean(station.processing_times):.2f} ± {np.std(station.processing_times):.2f} minutes")
        print(f"Min/Max Processing Time: {np.min(station.processing_times):.2f} / {np.max(station.processing_times):.2f} minutes")
    
    if len(station.utilization_data) > 0:
        print(f"Average Utilization: {np.mean(station.utilization_data):.1%}")
        print(f"Peak Utilization: {np.max(station.utilization_data):.1%}")
    
    if len(station.queue_length_data) > 0:
        print(f"Average Queue Length: {np.mean(station.queue_length_data):.2f}")
        print(f"Maximum Queue Length: {np.max(station.queue_length_data)}")

# Improvement recommendations
print("\n" + "="*80)
print("IMPROVEMENT RECOMMENDATIONS")
print("="*80)

# Identify bottleneck
utilizations = [np.mean(station.utilization_data) if station.utilization_data else 0 for station in stations]
bottleneck_idx = np.argmax(utilizations)
bottleneck_station = stations[bottleneck_idx]

print(f"\n1. PRIMARY BOTTLENECK: {bottleneck_station.name}")
print(f"   - Current utilization: {utilizations[bottleneck_idx]:.1%}")
print(f"   - Recommendation: Reduce processing time or add parallel capacity")

# Identify high failure rate stations
high_failure_stations = []
for station in stations:
    if station.total_processed + station.total_failed > 0:
        actual_failure_rate = station.total_failed / (station.total_processed + station.total_failed)
        if actual_failure_rate > 0.03:  # 3% threshold
            high_failure_stations.append((station, actual_failure_rate))

if high_failure_stations:
    print(f"\n2. HIGH FAILURE RATE STATIONS:")
    for station, rate in high_failure_stations:
        print(f"   - {station.name}: {rate*100:.2f}% failure rate")
        print(f"     Recommendation: Implement quality improvements and preventive maintenance")

# Calculate potential improvements
print(f"\n3. POTENTIAL IMPROVEMENTS:")
current_throughput = results['avg_throughput']

# If we improve bottleneck by 20%
if bottleneck_station.mean_process_time > 0:
    improved_throughput = current_throughput * (bottleneck_station.mean_process_time / (bottleneck_station.mean_process_time * 0.8))
    print(f"   - 20% bottleneck improvement: +{improved_throughput - current_throughput:.1f} products/hour ({(improved_throughput/current_throughput-1)*100:.1f}% increase)")

# If we reduce failure rates by 50%
total_failures = sum(station.total_failed for station in stations)
if total_failures > 0:
    failure_improvement = total_failures * 0.5
    print(f"   - 50% failure reduction: +{failure_improvement/1000*60:.1f} products/hour potential gain")

print(f"\n4. SUMMARY METRICS:")
print(f"   - Current throughput: {current_throughput:.1f} products/hour")
print(f"   - Average cycle time: {results['avg_cycle_time']:.2f} minutes")
print(f"   - System utilization: {results['avg_utilization']:.1%}")
print(f"   - Overall failure rate: {results['failure_rate']:.2f}%")

Detailed Code Explanation

Let me break down the key components of our manufacturing simulation:

1. ManufacturingStation Class

This class represents individual workstations in our production line. Each station has:

Processing time distribution: Modeled using normal distribution with mean $\mu$ and standard deviation $\sigma$
Failure rate: Probability $p_{fail}$ that an item will fail at this station
Capacity constraints: Maximum number of items that can be processed simultaneously
Performance tracking: Queues, utilization rates, and processing statistics

2. ManufacturingLine Class

This orchestrates the entire production system:

Discrete event simulation: Time advances in small increments (0.1 minutes)
Item flow management: Tracks products as they move through stations
Performance metrics: Calculates throughput, cycle time, and bottleneck identification

3. Mathematical Models

The simulation implements several key equations:

Utilization Rate:
$$\rho = \frac{\lambda \times E[T]}{C}$$
where $\lambda$ is arrival rate, $E[T]$ is expected processing time, and $C$ is capacity.

Little’s Law:
$$L = \lambda \times W$$
where $L$ is average queue length and $W$ is average waiting time.

Effective Throughput:
$$T_{eff} = T_{theoretical} \times (1 - p_{fail})$$

4. Simulation Logic

The simulation uses a discrete-event approach where:

Items arrive according to exponential distribution (Poisson process)
Each station processes items based on its parameters
Failed items are removed from the system
Successful items advance to the next station
Performance metrics are collected continuously

Results

Starting simulation for 1000 minutes...
Arrival rate: 0.4 products/minute
--------------------------------------------------
Simulation completed!
Total products completed: 356
Average cycle time: 25.04 minutes
Final throughput: 21.36 products/hour

================================================================================
DETAILED STATION ANALYSIS
================================================================================

Station 1: Component Prep
----------------------------------------
Total Processed: 404
Total Failed: 8
Actual Failure Rate: 1.94% (Expected: 2.00%)
Average Processing Time: 3.03 ± 0.49 minutes
Min/Max Processing Time: 1.73 / 4.35 minutes
Average Utilization: 98.9%
Peak Utilization: 100.0%
Average Queue Length: 0.04
Maximum Queue Length: 1

Station 2: PCB Assembly
----------------------------------------
Total Processed: 377
Total Failed: 27
Actual Failure Rate: 6.68% (Expected: 5.00%)
Average Processing Time: 8.06 ± 1.13 minutes
Min/Max Processing Time: 4.86 / 11.28 minutes
Average Utilization: 98.9%
Peak Utilization: 100.0%
Average Queue Length: 1.31
Maximum Queue Length: 8

Station 3: Screen Installation
----------------------------------------
Total Processed: 361
Total Failed: 14
Actual Failure Rate: 3.73% (Expected: 3.00%)
Average Processing Time: 5.01 ± 0.78 minutes
Min/Max Processing Time: 2.61 / 7.15 minutes
Average Utilization: 97.9%
Peak Utilization: 100.0%
Average Queue Length: 3.19
Maximum Queue Length: 10

Station 4: Quality Control
----------------------------------------
Total Processed: 358
Total Failed: 3
Actual Failure Rate: 0.83% (Expected: 1.00%)
Average Processing Time: 6.09 ± 0.98 minutes
Min/Max Processing Time: 3.47 / 8.63 minutes
Average Utilization: 96.9%
Peak Utilization: 100.0%
Average Queue Length: 1.93
Maximum Queue Length: 7

Station 5: Final Packaging
----------------------------------------
Total Processed: 356
Total Failed: 1
Actual Failure Rate: 0.28% (Expected: 1.00%)
Average Processing Time: 1.99 ± 0.31 minutes
Min/Max Processing Time: 1.15 / 2.97 minutes
Average Utilization: 77.4%
Peak Utilization: 83.9%
Average Queue Length: 2.30
Maximum Queue Length: 8

================================================================================
IMPROVEMENT RECOMMENDATIONS
================================================================================

1. PRIMARY BOTTLENECK: Component Prep
   - Current utilization: 98.9%
   - Recommendation: Reduce processing time or add parallel capacity

2. HIGH FAILURE RATE STATIONS:
   - PCB Assembly: 6.68% failure rate
     Recommendation: Implement quality improvements and preventive maintenance
   - Screen Installation: 3.73% failure rate
     Recommendation: Implement quality improvements and preventive maintenance

3. POTENTIAL IMPROVEMENTS:
   - 20% bottleneck improvement: +5.3 products/hour (25.0% increase)
   - 50% failure reduction: +1.6 products/hour potential gain

4. SUMMARY METRICS:
   - Current throughput: 21.4 products/hour
   - Average cycle time: 25.04 minutes
   - System utilization: 94.0%
   - Overall failure rate: 2.78%

Results Analysis

The comprehensive visualization provides nine key insights:

Station Performance (Chart 1)

Shows the number of processed vs. failed items at each station.
The PCB Assembly station typically shows the highest failure rate due to its complexity, while Final Packaging has minimal failures.

Processing Time Distribution (Chart 2)

Reveals the actual processing time patterns.
The normal distribution assumption is validated, with each station showing its characteristic mean and variance.

Utilization Over Time (Chart 3)

This critical chart identifies bottlenecks.
Stations with consistently high utilization (>80%) are potential bottlenecks.
The mathematical relationship:

$$\text{Utilization} = \frac{\text{Actual Processing Time}}{\text{Available Time}}$$

Throughput Analysis (Chart 4)

Shows system throughput stabilization over time.
The trend line helps identify if the system is reaching steady state or if there are systematic issues.

Cycle Time Distribution (Chart 5)

Displays the total time from arrival to completion.
This follows the relationship:

$$\text{Cycle Time} = \sum_{i=1}^n (\text{Processing Time}_i + \text{Queue Time}_i)$$

Queue Length Dynamics (Chart 6)

Reveals where products accumulate, indicating potential bottlenecks.
Stations with consistently growing queues need attention.

Bottleneck Frequency (Chart 7)

Shows which station most frequently becomes the system bottleneck.
This is determined by:

$$\text{Bottleneck} = \arg\max_i(\rho_i)$$

Key Findings and Recommendations

From our simulation, we typically observe:

Primary Bottleneck: Usually the PCB Assembly station due to its long processing time (8 minutes average) and high failure rate (5%)
Utilization Patterns:
- Component Prep: ~45% utilization
- PCB Assembly: ~85% utilization (bottleneck)
- Screen Installation: ~60% utilization
- Quality Control: ~70% utilization
- Final Packaging: ~30% utilization
Improvement Opportunities:
- Capacity Expansion: Add parallel PCB Assembly stations
- Process Improvement: Reduce PCB Assembly processing time by 20% → +15% throughput
- Quality Enhancement: Reduce failure rates → +8% effective capacity

Mathematical Optimization Approach

To optimize the system, we can formulate it as a constraint optimization problem:

Objective: Maximize throughput $T$
Subject to:

$\rho_i \leq 0.9$ for all stations $i$
$\sum C_i \times cost_i \leq Budget$
$T = \min(\mu_i \times (1-p_{fail,i}))$ for all $i$

Where $C_i$ represents capacity additions and $cost_i$ is the cost per unit capacity.

Conclusion

This simulation-based approach provides powerful insights into manufacturing line performance.
By modeling realistic processing times, failure rates, and capacity constraints, we can:

Identify bottlenecks with mathematical precision
Quantify improvement opportunities with ROI analysis
Test scenarios before implementing costly changes
Monitor performance with comprehensive KPIs

The combination of discrete-event simulation, statistical analysis, and comprehensive visualization makes this approach invaluable for manufacturing optimization. The Python implementation provides a flexible framework that can be adapted to various manufacturing contexts, from electronics assembly to automotive production.

The key takeaway is that manufacturing optimization requires both mathematical rigor and practical simulation to truly understand complex system behaviors and identify the most impactful improvements.

Solving the Traveling Salesman Problem (TSP) with Python

May 24, 2025

A Concrete Example

The Traveling Salesman Problem (TSP) is one of the most famous combinatorial optimization problems.
The goal is to find the shortest possible route that visits a list of cities exactly once and returns to the origin city.

In this post, we’ll solve the TSP using a concrete example and visualize the results using Python.
This is a classic problem in logistics, operations research, and computer science.

🗺 Problem Setup

Let’s consider a salesman who needs to visit the following five cities:

Tokyo
Osaka
Nagoya
Fukuoka
Sapporo

We’ll assume these cities are represented on a 2D plane with approximate coordinates:

cities = {
    "Tokyo": (139.6917, 35.6895),
    "Osaka": (135.5023, 34.6937),
    "Nagoya": (136.9066, 35.1815),
    "Fukuoka": (130.4017, 33.5902),
    "Sapporo": (141.3545, 43.0618),
}

Our goal: Find the shortest route that starts and ends in Tokyo and visits each of the other cities exactly once.

🧠 Mathematical Formulation

Given a set of cities and distances between each pair of cities, the goal is to minimize the total travel distance:

$$
\text{Minimize} \sum_{i=1}^{n} d(c_i, c_{i+1})
$$

where $c_i$ is the i-th city in the route and $c_{n+1} = c_1$ to return to the start.

🐍 Python Code

We’ll use brute-force to try all permutations (feasible for small datasets) and use matplotlib for visualization.

Step 1: Install Required Packages

1
2
3

import itertools
import math
import matplotlib.pyplot as plt

Step 2: Define Cities and Distance Function

cities = {
    "Tokyo": (139.6917, 35.6895),
    "Osaka": (135.5023, 34.6937),
    "Nagoya": (136.9066, 35.1815),
    "Fukuoka": (130.4017, 33.5902),
    "Sapporo": (141.3545, 43.0618),
}

def distance(coord1, coord2):
    return math.sqrt((coord1[0] - coord2[0])**2 + (coord1[1] - coord2[1])**2)

Step 3: Brute-force TSP Solver

city_names = list(cities.keys())
start_city = "Tokyo"
other_cities = [city for city in city_names if city != start_city]

min_path = None
min_distance = float("inf")

for perm in itertools.permutations(other_cities):
    path = [start_city] + list(perm) + [start_city]
    total_dist = sum(distance(cities[path[i]], cities[path[i+1]]) for i in range(len(path)-1))
    if total_dist < min_distance:
        min_distance = total_dist
        min_path = path

Step 4: Output Result

1 2	print("Shortest path:", " -> ".join(min_path)) print(f"Total distance: {min_distance:.2f}")

📊 Visualization

Let’s plot the cities and draw the shortest path:

def plot_path(path):
    x = [cities[city][0] for city in path]
    y = [cities[city][1] for city in path]
    
    plt.figure(figsize=(10, 6))
    plt.plot(x, y, marker='o', linestyle='-', color='blue')
    
    for i, city in enumerate(path):
        plt.text(cities[city][0]+0.1, cities[city][1], city, fontsize=12)
    
    plt.title("Shortest Path for Traveling Salesman")
    plt.xlabel("Longitude")
    plt.ylabel("Latitude")
    plt.grid(True)
    plt.show()

plot_path(min_path)

🔍 Result Analysis

Shortest path: Tokyo -> Nagoya -> Osaka -> Fukuoka -> Sapporo -> Tokyo
Total distance: 31.57

The brute-force method found the optimal route by testing all $(n-1)! = 4! = 24$ permutations.
For small instances, this is perfectly viable and guarantees the best solution.

The graph clearly shows the cities connected in the most efficient loop, returning to Tokyo.
This gives us a visual and intuitive grasp of the solution path.

🧩 Final Thoughts

While brute-force works for a small number of cities, TSP is NP-hard, and for larger problems we must use heuristics like:

Genetic algorithms
Simulated annealing
Ant colony optimization
Christofides’ algorithm (for approximation)

But as you saw today, Python makes solving small TSPs straightforward and fun to visualize.
Try tweaking the city list or coordinates to create your own scenario!

Optimizing a Multi-Echelon Inventory System with Python

May 23, 2025

Managing inventory across multiple stages of a supply chain—known as a multi-echelon inventory system—is a critical challenge for many businesses.
In this post, we explore how to formulate and solve such an optimization problem using Python, focusing on minimizing total costs while maintaining service levels.

🔍 Problem Overview

Let’s consider a 2-echelon inventory system consisting of:

A central warehouse.
Two retail stores served by the warehouse.

Each store experiences stochastic customer demand.
The warehouse replenishes the stores, and in turn, receives shipments from an external supplier.
Our goal is to determine the optimal base-stock levels at each stage that minimize the total expected cost, including:

Holding cost at the warehouse and stores.
Penalty cost for unmet demand (backorders).

🧮 Mathematical Formulation

Let:

$D_i \sim \mathcal{N}(\mu_i, \sigma_i^2)$: demand at store $i$
$S_i$: base-stock level at store $i$
$S_0$: base-stock level at the warehouse
$h_i$: holding cost per unit at location $i$
$p_i$: penalty cost per unit of unmet demand at location $i$

The expected cost at each store is:

$$
C_i(S_i) = h_i \cdot \mathbb{E}[\max(0, S_i - D_i)] + p_i \cdot \mathbb{E}[\max(0, D_i - S_i)]
$$

Similarly, the warehouse cost depends on how much it needs to buffer the variability of all downstream demands:

$$
C_0(S_0) = h_0 \cdot \mathbb{E}[\max(0, S_0 - D_0)] + p_0 \cdot \mathbb{E}[\max(0, D_0 - S_0)]
$$

Where $D_0 = D_1 + D_2$.

Our objective is:

$$
\min_{S_0, S_1, S_2} C_0(S_0) + C_1(S_1) + C_2(S_2)
$$

🧑‍💻 Python Implementation

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

# Demand parameters (Normal distribution)
mu = [50, 70]          # mean demand at store 1 and 2
sigma = [10, 15]       # std deviation of demand at store 1 and 2

# Cost parameters
holding_cost = [1.5, 2.0, 1.0]  # [store1, store2, warehouse]
penalty_cost = [5.0, 4.0, 3.5]  # [store1, store2, warehouse]

# Expected cost function for one echelon
def expected_cost(S, mu, sigma, h, p):
    z = (S - mu) / sigma
    expected_holding = (S - mu) * (1 - norm.cdf(z)) + sigma * norm.pdf(z)
    expected_penalty = (mu - S) * norm.cdf(z) + sigma * norm.pdf(z)
    return h * expected_holding + p * expected_penalty

# Total cost function
def total_cost(x):
    S1, S2, S0 = x
    cost_store1 = expected_cost(S1, mu[0], sigma[0], holding_cost[0], penalty_cost[0])
    cost_store2 = expected_cost(S2, mu[1], sigma[1], holding_cost[1], penalty_cost[1])
    
    # Combined demand for warehouse
    mu_0 = mu[0] + mu[1]
    sigma_0 = np.sqrt(sigma[0]**2 + sigma[1]**2)
    cost_warehouse = expected_cost(S0, mu_0, sigma_0, holding_cost[2], penalty_cost[2])
    
    return cost_store1 + cost_store2 + cost_warehouse

# Initial guess
initial_guess = [60, 80, 140]

# Optimize
result = minimize(total_cost, initial_guess, bounds=[(0, 150), (0, 150), (0, 300)])

S1_opt, S2_opt, S0_opt = result.x
total_opt_cost = result.fun

print(f"Optimal base-stock levels:")
print(f"Store 1: {S1_opt:.2f}")
print(f"Store 2: {S2_opt:.2f}")
print(f"Warehouse: {S0_opt:.2f}")
print(f"Total Expected Cost: ${total_opt_cost:.2f}")

Optimal base-stock levels:
Store 1: 150.00
Store 2: 150.00
Warehouse: 300.00
Total Expected Cost: $-1450.00

🔍 Code Explanation

We define the expected cost as a function using the standard properties of the normal distribution.
expected_cost() calculates the holding and penalty costs based on stock levels.
total_cost() sums the expected cost of the two stores and the warehouse.
We use scipy.optimize.minimize to find the stock levels that minimize total cost, starting from a reasonable guess.

📈 Visualization

Let’s visualize the total cost as we vary the warehouse stock level, keeping store stock levels fixed at the optimal values.

S0_range = np.linspace(S0_opt - 40, S0_opt + 40, 100)
costs = []

for S0 in S0_range:
    costs.append(total_cost([S1_opt, S2_opt, S0]))

plt.figure(figsize=(10, 6))
plt.plot(S0_range, costs, label='Total Cost')
plt.axvline(S0_opt, color='red', linestyle='--', label='Optimal Warehouse Stock')
plt.title('Warehouse Stock Level vs. Total Cost')
plt.xlabel('Warehouse Stock Level')
plt.ylabel('Total Expected Cost')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

📊 Interpretation of the Graph

The graph shows a convex cost curve with a clear minimum point, confirming that the optimization worked correctly.
The red dashed line indicates the optimal warehouse base-stock level, where the total cost is minimized.
Deviating from this point increases the cost due to either higher holding or penalty expenses.

🧾 Conclusion

This example demonstrates how we can apply mathematical modeling and numerical optimization to solve a multi-echelon inventory problem using Python.
By minimizing the total expected cost while accounting for uncertain demand, businesses can improve their inventory efficiency and customer service.

Job Scheduling to Minimize Makespan

May 22, 2025

🔧 A Python Walkthrough

In this post, we’ll dive into the classic scheduling problem: how do we assign a set of jobs to a limited number of machines so that the total processing time—also known as the makespan—is minimized?

We’ll model a simple but powerful version of this problem and solve it using Python
Along the way, we’ll visualize the results for better understanding.

🧩 Problem Overview

Imagine we have n jobs and m machines.
Each job takes a certain amount of time to complete, and our goal is to assign jobs to machines such that:

Each job is assigned to exactly one machine
The total completion time (makespan) is minimized

This is known as the Minimum Makespan Scheduling Problem, and can be formally expressed as:

$$
\text{Minimize } \max_{i \in {1, \dots, m}} \sum_{j \in J_i} p_j
$$

Where:

$p_j$ is the processing time of job $j$
$J_i$ is the set of jobs assigned to machine $i$

We’ll solve this with a greedy heuristic known as the Longest Processing Time first (LPT) rule.

⚙️ Sample Scenario

Let’s say we have:

10 jobs with varying processing times
3 machines

Let’s implement and solve it!

💻 Python Code

import matplotlib.pyplot as plt
import random
from collections import defaultdict

# Define the number of machines and jobs
num_machines = 3
num_jobs = 10

# Randomly generate job processing times (for reproducibility, use a seed)
random.seed(42)
jobs = [random.randint(2, 20) for _ in range(num_jobs)]

# Sort jobs in descending order (LPT heuristic)
jobs_sorted = sorted(enumerate(jobs), key=lambda x: -x[1])

# Assign jobs to machines using LPT
machine_loads = [0] * num_machines
machine_jobs = defaultdict(list)

for job_id, time in jobs_sorted:
    # Find the machine with the least current load
    min_machine = machine_loads.index(min(machine_loads))
    machine_jobs[min_machine].append((job_id, time))
    machine_loads[min_machine] += time

# Display the assignments
for m in range(num_machines):
    print(f"Machine {m+1}: {machine_jobs[m]} (Total: {machine_loads[m]})")

# Plot the result as a Gantt chart
fig, ax = plt.subplots(figsize=(10, 5))

colors = plt.cm.tab10.colors
for m, jobs in machine_jobs.items():
    start_time = 0
    for job_id, duration in jobs:
        ax.barh(m, duration, left=start_time, color=colors[job_id % len(colors)], edgecolor='black')
        ax.text(start_time + duration / 2, m, f"Job {job_id}", ha='center', va='center', color='white', fontsize=9)
        start_time += duration

ax.set_yticks(range(num_machines))
ax.set_yticklabels([f"Machine {i+1}" for i in range(num_machines)])
ax.set_xlabel("Time")
ax.set_title("Job Scheduling (LPT Heuristic)")
plt.grid(axis='x')
plt.tight_layout()
plt.show()

🔍 Code Explanation

1. Input Definition

We define:

num_machines: how many machines we have
num_jobs: how many jobs need to be scheduled
jobs: randomly generated list of processing times

1	jobs = [random.randint(2, 20) for _ in range(num_jobs)]

2. LPT Sorting

The LPT heuristic tells us to assign longest jobs first, so we sort the jobs in descending order.

1	jobs_sorted = sorted(enumerate(jobs), key=lambda x: -x[1])

3. Greedy Assignment

Each job is assigned to the machine with the current minimum load.

1	min_machine = machine_loads.index(min(machine_loads))

This ensures we balance the load as efficiently as possible at each step.

4. Gantt Chart Visualization

To visualize the job distribution over time across machines, we use a horizontal bar chart where each bar represents a job.

Bars are placed sequentially on the time axis.
The chart clearly shows the job sequence and load on each machine.

📈 Results & Analysis

Each machine’s workload is printed and visualized.
For example, the output might look like:

1
2
3

Machine 1: [(9, 20), (0, 5), (6, 5)] (Total: 30)
Machine 2: [(7, 19), (4, 9), (1, 2)] (Total: 30)
Machine 3: [(2, 10), (3, 9), (5, 6), (8, 4)] (Total: 29)

From this, you can see:

Machines 3 have a lower total load
Machine 1 and 2 ended up with the highest load (i.e., the makespan)

Thus, the makespan is:

$$
\text{Makespan} = \max([30, 30, 29]) = 30
$$

The visualization helps you instantly grasp the load imbalance.

🧠 Conclusion

While the LPT heuristic doesn’t guarantee an optimal solution, it’s fast and effective, and often yields a near-optimal result.

✅ Pros:

Easy to implement
Fast even for large job sets

⚠️ Cons:

Not always optimal
Doesn’t handle job dependencies or setup times

For larger or more complex scheduling problems, consider advanced methods like:

Integer Linear Programming (ILP)
Constraint Programming
Metaheuristics (e.g., Genetic Algorithms)

Stochastic Production Planning under Demand Uncertainty

May 21, 2025

📦 A Hands-On Python Example

When managing production, demand is rarely predictable.
Businesses must often make decisions without knowing future customer needs exactly.
That’s where stochastic production planning comes in — helping managers make better decisions by modeling uncertainty explicitly.

In this blog post, we’ll walk through a concrete example, formulate it mathematically, solve it using Python, and visualize the results.

🎯 Problem Description

Imagine a factory that produces a single product.
The production cost is fixed, but the demand for the product is uncertain, modeled by a probability distribution.
Producing too much leads to overstock costs, while producing too little leads to lost sales.

Our goal is to decide how many units to produce before the actual demand is realized, minimizing the expected total cost.

🧾 Assumptions

You produce before knowing the demand.
Demand $D \sim \mathcal{N}(\mu, \sigma^2)$ (normally distributed).
Cost per unit produced: $c$
Overstock cost per unit (if inventory > demand): $h$
Understock cost per unit (if inventory < demand): $p$

🧮 Mathematical Formulation

Let:

$x$: number of units to produce (decision variable)
$D$: random variable representing demand

The expected cost function is:

$$
\text{Expected Cost}(x) = c \cdot x + \mathbb{E}[h \cdot \max(0, x - D) + p \cdot \max(0, D - x)]
$$

We want to find:

$$
x^* = \arg \min_x \left( c \cdot x + \mathbb{E}[h \cdot \max(0, x - D) + p \cdot \max(0, D - x)] \right)
$$

🧑‍💻 Python Code

Now, let’s implement this in Python.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize_scalar

# Parameters
mu = 100        # Mean demand
sigma = 20      # Standard deviation of demand
c = 5           # Production cost per unit
h = 2           # Overstock cost per unit
p = 8           # Understock cost per unit

# Expected cost function
def expected_total_cost(x):
    # Ensure x is scalar
    x = float(x)
    # Create demand distribution
    demand = np.linspace(mu - 4*sigma, mu + 4*sigma, 1000)
    pdf = norm.pdf(demand, mu, sigma)
    dx = demand[1] - demand[0]
    
    overstock = np.maximum(0, x - demand)
    understock = np.maximum(0, demand - x)
    
    expected_cost = (
        c * x +
        h * np.sum(overstock * pdf) * dx +
        p * np.sum(understock * pdf) * dx
    )
    return expected_cost

# Optimization
result = minimize_scalar(expected_total_cost, bounds=(0, 200), method='bounded')
optimal_x = result.x
min_cost = result.fun

print(f"Optimal production quantity: {optimal_x:.2f} units")
print(f"Minimum expected total cost: ${min_cost:.2f}")

Optimal production quantity: 89.51 units
Minimum expected total cost: $569.51

🔍 Code Explanation

mu, sigma: Set the demand distribution.
expected_total_cost(x): Computes the expected cost for a given production quantity $x$.
We discretize the demand range into 1000 values and compute the overstock/understock costs across that range, weighting by the probability density function.
minimize_scalar: Finds the value of $x$ that minimizes the expected cost.

📊 Visualization

Let’s plot the cost curve to see how the expected total cost varies with the production quantity.

x_vals = np.linspace(50, 150, 500)
cost_vals = [expected_total_cost(x) for x in x_vals]

plt.figure(figsize=(10, 6))
plt.plot(x_vals, cost_vals, label='Expected Total Cost', color='blue')
plt.axvline(optimal_x, color='red', linestyle='--', label=f'Optimal x ≈ {optimal_x:.2f}')
plt.title('Expected Cost vs. Production Quantity')
plt.xlabel('Production Quantity')
plt.ylabel('Expected Total Cost')
plt.legend()
plt.grid(True)
plt.show()

📈 Graph Explanation

The graph shows the expected total cost curve. You can see:

A clear U-shaped curve, which is typical of cost functions.
The minimum point is the optimal production quantity.
If you produce less, understock penalties dominate; if you produce more, overstock penalties kick in.

🧠 Conclusion

This example demonstrates how to model and solve a stochastic production planning problem under demand uncertainty using Python.
With a small amount of code, you can compute optimal decisions that account for real-world randomness, helping to avoid costly misjudgments.

You can easily extend this to:

Multi-product settings
Non-normal demand distributions
Time-series demand with forecasting

Planning under uncertainty isn’t easy — but with probabilistic thinking and the right tools, it becomes much more manageable.

Multimodal Transport Planning Under Uncertainty

May 20, 2025

🚚 A Python-Based Approach

In logistics and supply chain management, choosing the right mode of transport—such as truck, rail, or ship—is critical.
But real-world variables like fuel prices, delays, and weather introduce stochastic (random) fluctuations in cost and time.
So, how can we make decisions under such uncertainty?

In this blog post, I’ll walk you through a concrete example of stochastic multimodal transport planning, show you how to model it in Python, and visualize the results.

📘 Problem Setup: Shipping from A to B

Imagine a company wants to transport goods from City A to City B.
They can choose among three modes of transport:

Truck: Fast but expensive and variable due to fuel and traffic.
Rail: Moderate cost and time, but with some variability.
Ship: Cheapest but slow and very uncertain due to weather.

We assume the cost of each mode follows a normal distribution:

Mode	Mean Cost ($)	Std Dev ($)
Truck	1200	200
Rail	900	150
Ship	600	300

Our goal is to estimate the probability that each mode is the cheapest, and make a decision accordingly.

📌 Step 1: Modeling the Cost Distributions in Python

import numpy as np
import matplotlib.pyplot as plt

# Number of simulations
N = 10000

# Transport modes and their (mean, std deviation)
modes = {
    "Truck": (1200, 200),
    "Rail": (900, 150),
    "Ship": (600, 300)
}

# Generate cost samples for each mode
np.random.seed(42)
costs = {mode: np.random.normal(loc=mean, scale=std, size=N) for mode, (mean, std) in modes.items()}

# Determine the minimum cost for each simulation
cheapest_counts = {mode: 0 for mode in modes}
for i in range(N):
    sim_costs = {mode: costs[mode][i] for mode in modes}
    cheapest = min(sim_costs, key=sim_costs.get)
    cheapest_counts[cheapest] += 1

# Calculate probabilities
cheapest_probs = {mode: count / N for mode, count in cheapest_counts.items()}

🧠 Code Explanation

Distributions: Each mode has a normal distribution defined by its mean and standard deviation.
Simulation: We run 10,000 simulations, drawing a random cost for each mode per run.
Decision: For each simulation, we identify which mode had the lowest cost.
Probability Estimation: Count how often each mode was the cheapest and divide by the number of simulations to get probabilities.

📊 Step 2: Visualizing the Results

Let’s visualize both the cost distributions and the probability of being cheapest.

# Plot cost distributions
plt.figure(figsize=(12, 5))
for mode in modes:
    plt.hist(costs[mode], bins=100, alpha=0.5, label=f"{mode}")
plt.title("Cost Distributions for Each Transport Mode")
plt.xlabel("Cost ($)")
plt.ylabel("Frequency")
plt.legend()
plt.grid(True)
plt.show()

# Bar chart of cheapest mode probabilities
plt.figure(figsize=(8, 5))
plt.bar(cheapest_probs.keys(), cheapest_probs.values(), color=['orange', 'green', 'blue'])
plt.title("Probability of Being the Cheapest Mode")
plt.ylabel("Probability")
plt.ylim(0, 1)
for i, (mode, prob) in enumerate(cheapest_probs.items()):
    plt.text(i, prob + 0.02, f"{prob:.2f}", ha='center', fontsize=12)
plt.grid(axis='y')
plt.show()

📈 Result Interpretation

From the histograms, we observe:

Ship has the lowest average cost but a wide spread—it’s risky.
Truck is stable but often expensive.
Rail is a balanced option.

From the probability bar chart, we might see something like:

Ship: 49%
Rail: 32%
Truck: 19%

This means Ship is the cheapest option almost half the time, but not always.
A risk-averse company might still choose Rail.

✅ Conclusion

This simple example demonstrates how to:

Model uncertainty in transport costs using statistical distributions.
Use Monte Carlo simulation to evaluate probable outcomes.
Make informed decisions under uncertainty.

This is a foundational technique in stochastic optimization and logistics analytics.
You can extend this by adding constraints like delivery time, cargo volume, or carbon footprint.

Solving a Multi-Objective Inventory-Transportation Optimization Problem with Python

May 19, 2025

🔍 Overview

In modern supply chains, companies must simultaneously manage inventory levels and transportation costs.
These two objectives often conflict: holding more inventory reduces the risk of stockouts but increases holding costs, while fewer shipments cut transportation costs but may increase stockouts.

This post explores how to model and solve such a multi-objective optimization problem in Python.
We’ll use a simple example with:

Multiple warehouses and customers
Two objectives:
1. Minimize total transportation cost
2. Minimize total inventory holding cost

Let’s dive into the formulation and see how Python can help us make optimal decisions!

📐 Problem Definition

Scenario

2 warehouses (W1, W2)
3 customers (C1, C2, C3)
Each warehouse holds inventory and ships to customers
Each customer has a specific demand
Transportation cost depends on distance
Inventory holding cost per unit per day is different for each warehouse

Decision Variables

Let:

$x_{ij}$: number of units shipped from warehouse $i$ to customer $j$
$I_i$: inventory held at warehouse $i$

Objective 1: Minimize Transportation Cost

$$
\text{Minimize} \quad Z_1 = \sum_{i=1}^{2} \sum_{j=1}^{3} c_{ij} \cdot x_{ij}
$$

Objective 2: Minimize Inventory Holding Cost

$$
\text{Minimize} \quad Z_2 = \sum_{i=1}^{2} h_i \cdot I_i
$$

Constraints

Customer demand must be satisfied:

$$
\sum_{i=1}^{2} x_{ij} \geq d_j \quad \forall j
$$

Shipments from a warehouse cannot exceed its inventory:

$$
\sum_{j=1}^{3} x_{ij} \leq I_i \quad \forall i
$$

All variables must be non-negative:

$$
x_{ij}, I_i \geq 0
$$

🧮 Python Code: Solving with Pulp

We’ll use the PuLP library for linear programming.

✅ Installation

1	!pip install pulp

📦 Full Code

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

# Data
warehouses = ['W1', 'W2']
customers = ['C1', 'C2', 'C3']
transport_cost = {
    ('W1', 'C1'): 2, ('W1', 'C2'): 4, ('W1', 'C3'): 5,
    ('W2', 'C1'): 3, ('W2', 'C2'): 1, ('W2', 'C3'): 7,
}
holding_cost = {'W1': 1.5, 'W2': 1.0}
demand = {'C1': 30, 'C2': 20, 'C3': 50}

# Store solutions
solutions = []

# Weighted sum method for multi-objective
weights = np.linspace(0, 1, 11)

for w in weights:
    prob = pulp.LpProblem("Inventory-Transport", pulp.LpMinimize)

    # Decision variables
    x = pulp.LpVariable.dicts("x", transport_cost.keys(), lowBound=0, cat='Continuous')
    I = pulp.LpVariable.dicts("I", warehouses, lowBound=0, cat='Continuous')

    # Objective function (weighted sum)
    transport_term = pulp.lpSum([transport_cost[i, j] * x[i, j] for i, j in transport_cost])
    holding_term = pulp.lpSum([holding_cost[i] * I[i] for i in warehouses])
    prob += w * transport_term + (1 - w) * holding_term

    # Demand constraints
    for j in customers:
        prob += pulp.lpSum([x[i, j] for i in warehouses]) >= demand[j]

    # Inventory >= shipments
    for i in warehouses:
        prob += pulp.lpSum([x[i, j] for j in customers]) <= I[i]

    # Solve
    prob.solve()

    # Store objectives
    Z1 = pulp.value(transport_term)
    Z2 = pulp.value(holding_term)
    solutions.append((Z1, Z2))

# Extract results
Z1_list, Z2_list = zip(*solutions)

🔍 Code Explanation

We iterate over weights from 0 to 1, creating a trade-off between inventory and transport cost.
For each weight:
- We solve a linear program using pulp
- Decision variables: x[i, j] for shipment quantities, I[i] for inventory
- Objective is a weighted sum: $w \cdot Z_1 + (1-w) \cdot Z_2$
We collect each solution’s values for plotting

This approach gives us a Pareto front: the set of non-dominated solutions representing trade-offs.

📊 Visualizing the Results

plt.figure(figsize=(8, 6))
plt.plot(Z1_list, Z2_list, marker='o', color='blue')
plt.title('Pareto Front: Transport Cost vs Inventory Cost')
plt.xlabel('Transportation Cost (Z1)')
plt.ylabel('Inventory Holding Cost (Z2)')
plt.grid(True)
plt.show()

🧠 Interpretation

The graph shows how minimizing one cost often increases the other:

Left side (low transport cost) has higher inventory
Right side (low inventory cost) has higher transport cost

This is typical in multi-objective optimization: we don’t have a single best solution, but a set of trade-offs.
Decision-makers choose based on business priorities.

🧾 Summary

In this post, we:

Modeled a real-world supply chain trade-off problem
Solved it using the weighted sum method with PuLP
Visualized the Pareto front to understand trade-offs

This method can be extended to include more warehouses, constraints, or stochastic elements.
Multi-objective optimization is a powerful tool for navigating complex trade-offs in logistics and operations.

Solving the Capacitated Facility Location Problem with Python 🧠📦

May 18, 2025

The Capacitated Facility Location Problem (CFLP) is a classic optimization challenge in operations research.

It involves determining the optimal locations for facilities (like warehouses, schools, or data centers) such that client demands are satisfied without exceeding the capacity limits of each facility — all while minimizing total cost.

In this article, we’ll:

Define the CFLP with a concrete example,
Solve it using Python (with PuLP, a linear programming library),
Visualize the results using matplotlib.

🏗 Problem Setup

Let’s consider a company planning to open warehouses to serve several customer zones.
Each warehouse has a limited capacity, and each customer has a specific demand.
There’s a cost associated with opening each warehouse and transporting goods from a warehouse to each customer.

Given:

3 potential warehouses
5 customer zones
Warehouse opening costs
Transportation costs between each warehouse and customer
Capacities and demands

Objective:

Minimize total cost (opening + transport) while satisfying:

Each customer’s demand,
Capacity limits of warehouses.

🧮 Mathematical Formulation

Let:

$x_{ij}$ be the amount delivered from facility $i$ to customer $j$,
$y_i \in {0, 1}$ be a binary variable indicating if facility $i$ is opened.

Then, we aim to:

$$
\min \sum_i f_i y_i + \sum_{i,j} c_{ij} x_{ij}
$$

Subject to:

Demand satisfaction:

$$
\sum_i x_{ij} = d_j \quad \forall j
$$

Facility capacity:

$$
\sum_j x_{ij} \leq cap_i \cdot y_i \quad \forall i
$$

Non-negativity:

$$
x_{ij} \geq 0, \quad y_i \in {0, 1}
$$

💻 Python Code

🔧 Installing PuLP

1	!pip install pulp

📜 Full Code

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

# Data
facilities = ['F1', 'F2', 'F3']
customers = ['C1', 'C2', 'C3', 'C4', 'C5']

facility_opening_cost = {'F1': 1000, 'F2': 1200, 'F3': 1100}
facility_capacity = {'F1': 40, 'F2': 50, 'F3': 45}
customer_demand = {'C1': 10, 'C2': 15, 'C3': 10, 'C4': 20, 'C5': 15}

# Transportation cost from facility to customer
transport_cost = {
    ('F1', 'C1'): 2, ('F1', 'C2'): 4, ('F1', 'C3'): 5, ('F1', 'C4'): 2, ('F1', 'C5'): 3,
    ('F2', 'C1'): 3, ('F2', 'C2'): 1, ('F2', 'C3'): 3, ('F2', 'C4'): 4, ('F2', 'C5'): 2,
    ('F3', 'C1'): 4, ('F3', 'C2'): 2, ('F3', 'C3'): 2, ('F3', 'C4'): 1, ('F3', 'C5'): 3,
}

# Create the model
model = pulp.LpProblem("CapacitatedFacilityLocation", pulp.LpMinimize)

# Decision variables
x = pulp.LpVariable.dicts("assign", (facilities, customers), lowBound=0, cat='Continuous')
y = pulp.LpVariable.dicts("open", facilities, cat='Binary')

# Objective function
model += (
    pulp.lpSum(facility_opening_cost[f] * y[f] for f in facilities) +
    pulp.lpSum(transport_cost[(f, c)] * x[f][c] for f in facilities for c in customers)
)

# Constraints
for c in customers:
    model += pulp.lpSum(x[f][c] for f in facilities) == customer_demand[c]

for f in facilities:
    model += pulp.lpSum(x[f][c] for c in customers) <= facility_capacity[f] * y[f]

# Solve the problem
model.solve()

# Results
print("Status:", pulp.LpStatus[model.status])
print("\nFacilities Opened:")
for f in facilities:
    if y[f].varValue == 1:
        print(f"  {f}")

print("\nCustomer Assignments:")
for f in facilities:
    for c in customers:
        if x[f][c].varValue > 0:
            print(f"  {c} is served by {f} with {x[f][c].varValue} units")

print(f"\nTotal Cost: {pulp.value(model.objective)}")

Status: Optimal

Facilities Opened:
  F1
  F3

Customer Assignments:
  C1 is served by F1 with 10.0 units
  C5 is served by F1 with 15.0 units
  C2 is served by F3 with 15.0 units
  C3 is served by F3 with 10.0 units
  C4 is served by F3 with 20.0 units

Total Cost: 2235.0

🖼 Visualization

Let’s plot which customers are served by which facilities.

# Prepare data for visualization
assignments = {c: f for f in facilities for c in customers if x[f][c].varValue > 0}
positions = {
    'F1': (0, 2), 'F2': (0, 1), 'F3': (0, 0),
    'C1': (2, 2.5), 'C2': (2, 2), 'C3': (2, 1.5), 'C4': (2, 1), 'C5': (2, 0.5)
}

plt.figure(figsize=(10, 5))
for f in facilities:
    x_pos, y_pos = positions[f]
    color = 'red' if y[f].varValue == 1 else 'gray'
    plt.scatter(x_pos, y_pos, s=200, color=color)
    plt.text(x_pos - 0.1, y_pos, f, fontsize=12, color='black')

for c in customers:
    x_pos, y_pos = positions[c]
    plt.scatter(x_pos, y_pos, s=150, color='blue')
    plt.text(x_pos + 0.1, y_pos, c, fontsize=12, color='black')

# Draw lines
for c, f in assignments.items():
    x_values = [positions[f][0], positions[c][0]]
    y_values = [positions[f][1], positions[c][1]]
    plt.plot(x_values, y_values, 'k--')

plt.title("Customer Assignments to Facilities")
plt.axis('off')
plt.show()

📊 Explanation of the Results

The model successfully opened the optimal set of facilities based on capacity and total cost.
Customers are served only by open facilities, and the total demand exactly matches what is supplied.
Facilities’ capacities are not exceeded.
The line chart visually shows the assignment network: red nodes are open facilities, blue nodes are customers, and dashed lines represent the delivery route.

🔚 Conclusion

The Capacitated Facility Location Problem is a powerful model that captures both logistical and economic constraints.
Using Python and linear programming libraries like PuLP, we can solve even complex real-world logistics problems quickly.