Real-World Structural Optimization with Python

June 6, 2025

🏗️ Minimize Beam Weight Under Stress Constraints

📌 Introduction

Structural optimization is a powerful tool in engineering that aims to design structures that are both efficient and safe. One practical and commonly encountered example is minimizing the weight of a structural beam while ensuring it can withstand a given load without yielding.

In this post, we’ll walk through a real example of optimizing the cross-sectional dimensions of a cantilever beam under a point load. We’ll do this using Python and visualize the results clearly.

🎯 Problem Statement

We aim to minimize the weight of a rectangular cantilever beam with length $L$, subjected to a point load $P$ at the free end. The beam is made of steel, and we must ensure that the maximum bending stress does not exceed the yield stress of the material.

Variables:

Width of cross-section: $b$
Height of cross-section: $h$ (design variable)

Given:

Length $L = 2 , \text{m}$
Load $P = 1000 , \text{N}$
Density of steel $\rho = 7850 , \text{kg/m}^3$
Yield stress $\sigma_{\text{yield}} = 250 \times 10^6 , \text{Pa}$

Objective:

Minimize the weight:

$$
\text{Weight} = \rho \cdot g \cdot (b \cdot h \cdot L)
$$

Subject to:

$$
\sigma_{\text{max}} = \frac{6PL}{b h^2} \leq \sigma_{\text{yield}}
$$

🔧 Python Code

Let’s solve this optimization problem using scipy.optimize.

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

# Given constants
P = 1000          # Load in N
L = 2.0           # Length in meters
rho = 7850        # Density of steel (kg/m^3)
g = 9.81          # Gravity (m/s^2)
sigma_yield = 250e6  # Yield stress in Pa
b = 0.05          # Fixed width in meters

# Objective function: weight of the beam
def weight(x):
    h = x[0]
    volume = b * h * L
    return rho * g * volume  # Weight = mass * gravity

# Constraint: maximum bending stress must be less than yield
def stress_constraint(x):
    h = x[0]
    sigma_max = 6 * P * L / (b * h**2)
    return sigma_yield - sigma_max

# Bounds and constraints
bounds = [(0.01, 0.5)]  # height h should be between 1cm and 50cm
constraints = [{'type': 'ineq', 'fun': stress_constraint}]

# Initial guess
x0 = [0.1]

# Run optimization
result = minimize(weight, x0, method='SLSQP', bounds=bounds, constraints=constraints)

# Extract result
h_opt = result.x[0]
w_opt = weight([h_opt])

print(f"Optimal height h = {h_opt:.4f} m")
print(f"Minimum weight = {w_opt:.2f} N")

🔍 Code Explanation

🧮 Objective Function

def weight(x):
    h = x[0]
    volume = b * h * L
    return rho * g * volume

We compute the volume of the beam and then convert it into weight by multiplying with density and gravity.

⚠️ Stress Constraint

def stress_constraint(x):
    h = x[0]
    sigma_max = 6 * P * L / (b * h**2)
    return sigma_yield - sigma_max

This ensures that the stress due to the load doesn’t exceed the material’s yield stress. The constraint returns positive values when valid.

🧠 Optimization

We use SLSQP, a sequential quadratic programming solver, with bounds and nonlinear constraints. We optimize the beam’s height while keeping the width fixed.

📈 Result Visualization

Let’s visualize how stress and weight vary with height.

h_vals = np.linspace(0.01, 0.5, 100)
weights = rho * g * b * h_vals * L
stresses = 6 * P * L / (b * h_vals**2)

fig, ax1 = plt.subplots()

ax1.set_xlabel('Height h (m)')
ax1.set_ylabel('Weight (N)', color='tab:blue')
ax1.plot(h_vals, weights, label='Weight', color='tab:blue')
ax1.tick_params(axis='y', labelcolor='tab:blue')

ax2 = ax1.twinx()
ax2.set_ylabel('Max Bending Stress (Pa)', color='tab:red')
ax2.plot(h_vals, stresses, label='Stress', color='tab:red')
ax2.axhline(sigma_yield, color='gray', linestyle='--', label='Yield Stress')
ax2.tick_params(axis='y', labelcolor='tab:red')

plt.title("Beam Weight and Stress vs Height")
fig.tight_layout()
plt.legend()
plt.grid(True)
plt.show()

Optimal height h = 0.0310 m
Minimum weight = 238.60 N

📊 Graph Interpretation

Blue curve: Shows how weight increases with height. Larger cross-sections are heavier.
Red curve: Shows how bending stress decreases with increasing height. Taller beams are stronger in bending.
The optimal solution is a balance where the stress just touches the yield limit — no material is wasted.

✅ Conclusion

This example illustrates a fundamental trade-off in structural engineering: minimizing weight while maintaining safety. Using Python and a few powerful libraries, we were able to:

Define a physics-based optimization problem
Solve it efficiently with real-world parameters
Visualize and interpret the solution

Protein Folding Optimization

June 5, 2025

A Practical Deep Dive with Python

Protein folding is one of the most fascinating challenges in computational biology. Today, we’ll explore a realistic protein folding optimization problem using Python, focusing on the HP (Hydrophobic-Polar) model - a simplified yet powerful approach to understanding protein structure.

The Problem: HP Model Protein Folding

The HP model represents proteins as sequences of hydrophobic (H) and polar (P) amino acids on a 2D lattice. The goal is to find the optimal folding configuration that minimizes energy by maximizing hydrophobic contacts while avoiding steric clashes.

Our objective function is:
$$E = -\sum_{i,j} \delta_{H_i H_j} \cdot \text{contact}(i,j)$$

Where $\delta_{H_i H_j} = 1$ if both residues $i$ and $j$ are hydrophobic and adjacent (but not consecutive in sequence), and 0 otherwise.

Let’s solve this step by step with a realistic example!

import numpy as np
import matplotlib.pyplot as plt
import random
from collections import defaultdict
import seaborn as sns
from scipy.optimize import differential_evolution
import warnings
warnings.filterwarnings('ignore')

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

class HPProteinFolder:
    """
    HP Model Protein Folding Optimizer
    
    This class implements a simplified protein folding model where:
    - H represents hydrophobic amino acids
    - P represents polar amino acids
    - Goal: Maximize hydrophobic contacts while avoiding collisions
    """
    
    def __init__(self, sequence):
        """
        Initialize the protein folder with an HP sequence
        
        Args:
            sequence (str): String of 'H' and 'P' characters representing the protein
        """
        self.sequence = sequence
        self.length = len(sequence)
        self.moves = [(0, 1), (1, 0), (0, -1), (-1, 0)]  # Up, Right, Down, Left
        self.move_names = ['U', 'R', 'D', 'L']
        
    def decode_folding(self, angles):
        """
        Convert angle representation to 2D coordinates
        
        Args:
            angles (list): List of angles (0-3) representing folding directions
            
        Returns:
            tuple: (coordinates, valid_folding_flag)
        """
        if len(angles) != self.length - 1:
            return None, False
            
        coordinates = [(0, 0)]  # Start at origin
        current_pos = (0, 0)
        occupied = {(0, 0)}
        
        for i, angle in enumerate(angles):
            # Convert continuous angle to discrete move
            move_idx = int(angle) % 4
            dx, dy = self.moves[move_idx]
            new_pos = (current_pos[0] + dx, current_pos[1] + dy)
            
            # Check for collision (self-intersection)
            if new_pos in occupied:
                return None, False
                
            coordinates.append(new_pos)
            occupied.add(new_pos)
            current_pos = new_pos
            
        return coordinates, True
    
    def calculate_energy(self, coordinates):
        """
        Calculate the energy of a folding configuration
        
        Energy = -1 * (number of H-H contacts that are adjacent in space but not in sequence)
        
        Args:
            coordinates (list): List of (x, y) coordinates for each amino acid
            
        Returns:
            float: Energy value (lower is better)
        """
        if coordinates is None:
            return float('inf')  # Invalid folding has infinite energy
            
        energy = 0
        
        # Check all pairs of amino acids
        for i in range(self.length):
            for j in range(i + 2, self.length):  # Skip adjacent in sequence
                # Only count H-H interactions
                if self.sequence[i] == 'H' and self.sequence[j] == 'H':
                    # Check if they are adjacent in 2D space
                    dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                    if dist == 1:  # Manhattan distance = 1 means adjacent
                        energy -= 1  # Favorable interaction
                        
        return energy
    
    def objective_function(self, angles):
        """
        Objective function for optimization
        
        Args:
            angles (array): Array of folding angles
            
        Returns:
            float: Energy to minimize
        """
        coordinates, valid = self.decode_folding(angles)
        if not valid:
            return 1000  # Heavy penalty for invalid folding
            
        return self.calculate_energy(coordinates)
    
    def optimize_folding(self, method='differential_evolution', max_iterations=1000):
        """
        Optimize protein folding using specified optimization method
        
        Args:
            method (str): Optimization method to use
            max_iterations (int): Maximum number of iterations
            
        Returns:
            tuple: (best_angles, best_energy, optimization_history)
        """
        bounds = [(0, 3.99) for _ in range(self.length - 1)]
        history = []
        
        def callback(xk, convergence=None):
            energy = self.objective_function(xk)
            history.append(energy)
            return False
        
        if method == 'differential_evolution':
            result = differential_evolution(
                self.objective_function,
                bounds,
                maxiter=max_iterations // 10,
                popsize=15,
                callback=callback,
                seed=42
            )
            best_angles = result.x
            best_energy = result.fun
            
        return best_angles, best_energy, history
    
    def monte_carlo_optimization(self, max_iterations=10000, temperature=1.0, cooling_rate=0.995):
        """
        Monte Carlo optimization with simulated annealing
        
        Args:
            max_iterations (int): Maximum number of iterations
            temperature (float): Initial temperature
            cooling_rate (float): Temperature cooling rate
            
        Returns:
            tuple: (best_angles, best_energy, history)
        """
        current_angles = [random.uniform(0, 3.99) for _ in range(self.length - 1)]
        current_energy = self.objective_function(current_angles)
        
        best_angles = current_angles.copy()
        best_energy = current_energy
        history = [current_energy]
        
        temp = temperature
        
        for iteration in range(max_iterations):
            # Generate neighbor solution
            new_angles = current_angles.copy()
            idx = random.randint(0, len(new_angles) - 1)
            new_angles[idx] = random.uniform(0, 3.99)
            
            new_energy = self.objective_function(new_angles)
            
            # Accept or reject based on Metropolis criterion
            if new_energy < current_energy or random.random() < np.exp(-(new_energy - current_energy) / temp):
                current_angles = new_angles
                current_energy = new_energy
                
                if current_energy < best_energy:
                    best_angles = current_angles.copy()
                    best_energy = current_energy
            
            history.append(best_energy)
            temp *= cooling_rate
            
            if iteration % 1000 == 0:
                print(f"Iteration {iteration}: Best Energy = {best_energy:.3f}, Temp = {temp:.4f}")
        
        return best_angles, best_energy, history

# Example protein sequences for testing
sequences = {
    'simple': 'HPHPPHHPHH',  # 10 amino acids
    'medium': 'HPHPPHHPHHPPHPPHHPHH',  # 20 amino acids  
    'complex': 'HHPPHPHPHPPHPPHPPHPHPHHPPHPPH'  # 27 amino acids
}

print("=== Protein Folding Optimization Results ===\n")

# Let's work with the medium complexity sequence
sequence = sequences['medium']
print(f"Protein sequence: {sequence}")
print(f"Length: {len(sequence)} amino acids")
print(f"Hydrophobic residues: {sequence.count('H')}")
print(f"Polar residues: {sequence.count('P')}\n")

# Initialize the protein folder
folder = HPProteinFolder(sequence)

# Optimize using differential evolution
print("Running Differential Evolution optimization...")
de_angles, de_energy, de_history = folder.optimize_folding(method='differential_evolution')
de_coords, de_valid = folder.decode_folding(de_angles)

print(f"DE - Best energy: {de_energy:.3f}")
print(f"DE - Valid folding: {de_valid}\n")

# Optimize using Monte Carlo with simulated annealing
print("Running Monte Carlo with Simulated Annealing...")
mc_angles, mc_energy, mc_history = folder.monte_carlo_optimization(max_iterations=5000)
mc_coords, mc_valid = folder.decode_folding(mc_angles)

print(f"MC - Best energy: {mc_energy:.3f}")
print(f"MC - Valid folding: {mc_valid}\n")

# Analysis and visualization
def analyze_folding(coordinates, sequence, title):
    """
    Analyze and display folding statistics
    """
    if coordinates is None:
        print(f"{title}: Invalid folding")
        return
        
    # Calculate statistics
    h_positions = [i for i, aa in enumerate(sequence) if aa == 'H']
    p_positions = [i for i, aa in enumerate(sequence) if aa == 'P']
    
    # Count H-H contacts
    hh_contacts = 0
    for i in range(len(sequence)):
        for j in range(i + 2, len(sequence)):
            if sequence[i] == 'H' and sequence[j] == 'H':
                dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                if dist == 1:
                    hh_contacts += 1
    
    # Calculate compactness (radius of gyration)
    center_x = np.mean([coord[0] for coord in coordinates])
    center_y = np.mean([coord[1] for coord in coordinates])
    rg = np.sqrt(np.mean([(coord[0] - center_x)**2 + (coord[1] - center_y)**2 for coord in coordinates]))
    
    print(f"{title} Analysis:")
    print(f"  H-H contacts: {hh_contacts}")
    print(f"  Energy: {-hh_contacts}")
    print(f"  Radius of gyration: {rg:.3f}")
    print(f"  Span X: {max(coord[0] for coord in coordinates) - min(coord[0] for coord in coordinates)}")
    print(f"  Span Y: {max(coord[1] for coord in coordinates) - min(coord[1] for coord in coordinates)}")
    print()

# Analyze both solutions
analyze_folding(de_coords, sequence, "Differential Evolution")
analyze_folding(mc_coords, sequence, "Monte Carlo")

# Create comprehensive visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Protein Folding Optimization Results', fontsize=16, fontweight='bold')

# Plot 1: DE Folding Structure
if de_coords and de_valid:
    ax = axes[0, 0]
    x_coords = [coord[0] for coord in de_coords]
    y_coords = [coord[1] for coord in de_coords]
    
    # Plot backbone
    ax.plot(x_coords, y_coords, 'k-', linewidth=2, alpha=0.7, label='Backbone')
    
    # Plot amino acids with different colors
    for i, (x, y) in enumerate(de_coords):
        color = 'red' if sequence[i] == 'H' else 'blue'
        marker = 'o' if sequence[i] == 'H' else 's'
        ax.scatter(x, y, c=color, s=150, marker=marker, 
                  edgecolors='black', linewidth=1, zorder=5)
        ax.annotate(f'{i}', (x, y), xytext=(5, 5), textcoords='offset points', 
                   fontsize=8, ha='left')
    
    ax.set_title(f'Differential Evolution\nEnergy: {de_energy:.1f}', fontweight='bold')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.grid(True, alpha=0.3)
    ax.legend(['Backbone', 'Hydrophobic (H)', 'Polar (P)'])
    ax.set_aspect('equal')

# Plot 2: MC Folding Structure  
if mc_coords and mc_valid:
    ax = axes[0, 1]
    x_coords = [coord[0] for coord in mc_coords]
    y_coords = [coord[1] for coord in mc_coords]
    
    # Plot backbone
    ax.plot(x_coords, y_coords, 'k-', linewidth=2, alpha=0.7)
    
    # Plot amino acids
    for i, (x, y) in enumerate(mc_coords):
        color = 'red' if sequence[i] == 'H' else 'blue'
        marker = 'o' if sequence[i] == 'H' else 's'
        ax.scatter(x, y, c=color, s=150, marker=marker, 
                  edgecolors='black', linewidth=1, zorder=5)
        ax.annotate(f'{i}', (x, y), xytext=(5, 5), textcoords='offset points', 
                   fontsize=8, ha='left')
    
    ax.set_title(f'Monte Carlo\nEnergy: {mc_energy:.1f}', fontweight='bold')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

# Plot 3: Convergence Comparison
ax = axes[0, 2]
if len(de_history) > 0:
    ax.plot(range(len(de_history)), de_history, 'g-', linewidth=2, label='Differential Evolution')
if len(mc_history) > 0:
    # Sample MC history for better visualization
    sample_indices = np.linspace(0, len(mc_history)-1, min(200, len(mc_history)), dtype=int)
    sampled_history = [mc_history[i] for i in sample_indices]
    ax.plot(sample_indices, sampled_history, 'r-', linewidth=2, label='Monte Carlo')

ax.set_title('Optimization Convergence', fontweight='bold')
ax.set_xlabel('Iteration')
ax.set_ylabel('Best Energy Found')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 4: Energy Landscape Analysis
ax = axes[1, 0]
# Generate random foldings to show energy distribution
random_energies = []
for _ in range(1000):
    random_angles = [random.uniform(0, 3.99) for _ in range(len(sequence) - 1)]
    energy = folder.objective_function(random_angles)
    if energy < 100:  # Only valid foldings
        random_energies.append(energy)

if random_energies:
    ax.hist(random_energies, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
    ax.axvline(de_energy, color='green', linestyle='--', linewidth=2, label=f'DE: {de_energy:.1f}')
    ax.axvline(mc_energy, color='red', linestyle='--', linewidth=2, label=f'MC: {mc_energy:.1f}')
    ax.set_title('Energy Distribution\n(Random vs Optimized)', fontweight='bold')
    ax.set_xlabel('Energy')
    ax.set_ylabel('Frequency')
    ax.legend()
    ax.grid(True, alpha=0.3)

# Plot 5: Sequence Analysis
ax = axes[1, 1]
positions = list(range(len(sequence)))
colors = ['red' if aa == 'H' else 'blue' for aa in sequence]
bars = ax.bar(positions, [1]*len(sequence), color=colors, alpha=0.7, edgecolor='black')

ax.set_title('Protein Sequence', fontweight='bold')
ax.set_xlabel('Position')
ax.set_ylabel('Amino Acid Type')
ax.set_xticks(positions[::2])
ax.set_xticklabels(positions[::2])
ax.set_yticks([0.5, 1])
ax.set_yticklabels(['', 'H/P'])

# Add legend
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='red', alpha=0.7, label='Hydrophobic (H)'),
                  Patch(facecolor='blue', alpha=0.7, label='Polar (P)')]
ax.legend(handles=legend_elements)

# Plot 6: Contact Map
ax = axes[1, 2]
contact_matrix = np.zeros((len(sequence), len(sequence)))

# Fill contact matrix for the best solution (MC in this case)
best_coords = mc_coords if mc_valid else de_coords
if best_coords:
    for i in range(len(sequence)):
        for j in range(len(sequence)):
            if i != j:
                dist = abs(best_coords[i][0] - best_coords[j][0]) + abs(best_coords[i][1] - best_coords[j][1])
                if dist == 1:  # Adjacent in space
                    contact_matrix[i, j] = 1

im = ax.imshow(contact_matrix, cmap='RdYlBu_r', aspect='equal')
ax.set_title('Contact Map\n(Best Solution)', fontweight='bold')
ax.set_xlabel('Residue Index')
ax.set_ylabel('Residue Index')

# Add colorbar
cbar = plt.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('Contact')

plt.tight_layout()
plt.show()

# Performance comparison
print("=== Performance Comparison ===")
print(f"Differential Evolution:")
print(f"  Final Energy: {de_energy:.3f}")
print(f"  Convergence Steps: {len(de_history)}")

print(f"Monte Carlo + Simulated Annealing:")
print(f"  Final Energy: {mc_energy:.3f}")
print(f"  Total Steps: {len(mc_history)}")

# Determine the winner
if de_valid and mc_valid:
    if de_energy < mc_energy:
        print(f"\n🏆 Winner: Differential Evolution (Energy: {de_energy:.3f})")
    elif mc_energy < de_energy:
        print(f"\n🏆 Winner: Monte Carlo (Energy: {mc_energy:.3f})")
    else:
        print(f"\n🤝 Tie! Both methods achieved energy: {de_energy:.3f}")
else:
    print(f"\n⚠️  Some optimizations produced invalid foldings")

print("\n=== Mathematical Analysis ===")
print("Energy Function:")
print("E = -∑(i,j) δ_HiHj * contact(i,j)")
print("where δ_HiHj = 1 if both residues i,j are hydrophobic and adjacent")
print("Lower energy indicates better folding (more H-H contacts)")

Code Breakdown and Explanation

Let me walk you through the key components of this protein folding optimization implementation:

1. HPProteinFolder Class Architecture

The core class implements the HP model with several critical methods:

__init__: Initializes the protein sequence and defines possible moves (Up, Right, Down, Left) on a 2D lattice
decode_folding: Converts optimization variables (angles) into 2D coordinates while checking for self-intersections
calculate_energy: Implements our energy function $E = -\sum_{i,j} \delta_{H_i H_j} \cdot \text{contact}(i,j)$

2. Energy Function Deep Dive

The energy calculation is the heart of our optimization:

def calculate_energy(self, coordinates):
    energy = 0
    for i in range(self.length):
        for j in range(i + 2, self.length):  # Skip adjacent in sequence
            if self.sequence[i] == 'H' and self.sequence[j] == 'H':
                dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                if dist == 1:  # Manhattan distance = 1 means adjacent
                    energy -= 1  # Favorable interaction
    return energy

This implements the mathematical model where we:

Only consider hydrophobic-hydrophobic (H-H) interactions
Skip consecutive amino acids in the sequence (i+2 constraint)
Use Manhattan distance to determine spatial adjacency
Assign -1 energy per H-H contact (lower energy = better folding)

3. Optimization Algorithms

We implement two complementary approaches:

Differential Evolution: A population-based global optimizer that:

Maintains a population of candidate solutions
Uses mutation, crossover, and selection operations
Excellent for finding global optima in continuous spaces

Monte Carlo with Simulated Annealing: A probabilistic method that:

Accepts worse solutions with probability $P = e^{-\Delta E/T}$
Gradually reduces temperature $T$ to focus the search
Balances exploration and exploitation

4. Constraint Handling

The decode_folding method ensures valid protein structures by:

Tracking occupied lattice positions
Rejecting self-intersecting conformations
Returning infinite energy for invalid foldings

Results Analysis and Interpretation

When you run this code, you’ll observe several key insights:

=== Protein Folding Optimization Results ===

Protein sequence: HPHPPHHPHHPPHPPHHPHH
Length: 20 amino acids
Hydrophobic residues: 11
Polar residues: 9

Running Differential Evolution optimization...
DE - Best energy: -6.000
DE - Valid folding: True

Running Monte Carlo with Simulated Annealing...
Iteration 0: Best Energy = 1000.000, Temp = 0.9950
Iteration 1000: Best Energy = -4.000, Temp = 0.0066
Iteration 2000: Best Energy = -4.000, Temp = 0.0000
Iteration 3000: Best Energy = -4.000, Temp = 0.0000
Iteration 4000: Best Energy = -4.000, Temp = 0.0000
MC - Best energy: -4.000
MC - Valid folding: True

Differential Evolution Analysis:
  H-H contacts: 6
  Energy: -6
  Radius of gyration: 2.599
  Span X: 2
  Span Y: 8

Monte Carlo Analysis:
  H-H contacts: 4
  Energy: -4
  Radius of gyration: 2.035
  Span X: 4
  Span Y: 5

=== Performance Comparison ===
Differential Evolution:
  Final Energy: -6.000
  Convergence Steps: 100
Monte Carlo + Simulated Annealing:
  Final Energy: -4.000
  Total Steps: 5001

🏆 Winner: Differential Evolution (Energy: -6.000)

=== Mathematical Analysis ===
Energy Function:
E = -∑(i,j) δ_HiHj * contact(i,j)
where δ_HiHj = 1 if both residues i,j are hydrophobic and adjacent
Lower energy indicates better folding (more H-H contacts)

Optimization Performance

The algorithms typically find folding configurations with energies between -2 and -6, depending on the sequence. The negative values indicate successful H-H contact formation.

Convergence Behavior

Differential Evolution shows smooth, monotonic convergence
Monte Carlo exhibits more volatile behavior initially, then stabilizes as temperature decreases

Structural Analysis

The visualization reveals:

Red circles: Hydrophobic residues (H) that drive folding
Blue squares: Polar residues (P) that prefer solvent exposure
Black lines: Protein backbone connectivity
Contact maps: Spatial adjacency patterns

Energy Landscape

The histogram shows that most random foldings have poor energy (around 0 to -2), while our optimized solutions achieve significantly better values (typically -4 to -6).

Mathematical Insights

The optimization problem can be formulated as:

$$\min_{\theta} E(\theta) = \min_{\theta} \left( -\sum_{i < j-1} H_i H_j \cdot \text{adjacent}(pos_i(\theta), pos_j(\theta)) \right)$$

Where:

$\theta$ represents folding angles
$H_i \in {0,1}$ indicates if residue $i$ is hydrophobic
$pos_i(\theta)$ gives the 2D coordinates of residue $i$
$\text{adjacent}(p_1, p_2) = 1$ if $|p_1 - p_2|_1 = 1$

This is a challenging combinatorial optimization problem because:

Constraint complexity: Valid foldings must avoid self-intersection
Discrete-continuous hybrid: Lattice positions are discrete, but we optimize over continuous angles
Multimodal landscape: Many local optima exist

Practical Applications

This HP model, while simplified, captures essential protein folding principles:

Hydrophobic collapse: Nonpolar residues cluster together
Geometric constraints: Physical impossibility of self-intersection
Energy minimization: Nature favors thermodynamically stable conformations

Real protein folding involves additional factors like hydrogen bonding, electrostatic interactions, and entropy, but the HP model provides valuable insights into the fundamental optimization challenge that evolution has solved through natural selection.

The techniques demonstrated here scale to larger proteins and can be extended with more sophisticated energy functions, making this a practical foundation for computational structural biology applications.

Traffic Signal Optimization with Python

June 4, 2025

🚦 A Realistic Example

Urban traffic congestion is a growing problem. One way to alleviate this issue is through intelligent traffic signal control. In this blog post, we will explore a simplified yet realistic traffic signal optimization problem using Python. We’ll focus on minimizing total vehicle waiting time at a 4-way intersection and visualize the result to understand how well our model performs.

🧩 Problem Overview: 4-Way Intersection

Imagine a simple 4-way intersection with traffic coming from North-South (NS) and East-West (EW). Each direction has a green light duration, and vehicles arrive at different rates.

Our goal is:

Minimize total average waiting time by adjusting green light durations.

Assumptions:

Vehicles arrive following a Poisson distribution.
Signal cycle time is fixed (e.g., 60 seconds).
Only one direction can be green at a time.
The green duration for NS is g, and for EW is 60 - g.

📐 Mathematical Model

Let:

$\lambda_{NS}$: arrival rate of vehicles from North-South per second.
$\lambda_{EW}$: arrival rate from East-West per second.
$g$: green light time for NS direction.
$C$: total cycle time, fixed at 60 seconds.

Each queue has a waiting time modeled by:

$$
W = \text{arrival rate} \times \text{red time}^2 / (2 \times \text{green time})
$$

Total cost function (to minimize):

$$
\text{Total Waiting Time} = \frac{\lambda_{NS} \cdot (C - g)^2}{2g} + \frac{\lambda_{EW} \cdot g^2}{2(C - g)}
$$

🧪 Let’s Solve It in Python

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

# Constants
CYCLE_TIME = 60  # Total signal cycle time (seconds)
lambda_NS = 0.6  # Arrival rate per second from North-South
lambda_EW = 0.4  # Arrival rate per second from East-West

# Cost function to minimize: total waiting time
def total_waiting_time(g):
    if g <= 0 or g >= CYCLE_TIME:
        return np.inf  # avoid invalid green times
    red_NS = CYCLE_TIME - g
    red_EW = g
    wait_NS = (lambda_NS * red_NS ** 2) / (2 * g)
    wait_EW = (lambda_EW * red_EW ** 2) / (2 * (CYCLE_TIME - g))
    return wait_NS + wait_EW

# Find the optimal green time for NS direction
result = minimize_scalar(total_waiting_time, bounds=(1, CYCLE_TIME - 1), method='bounded')

optimal_g = result.x
min_wait = result.fun

print(f"✅ Optimal green time for NS: {optimal_g:.2f} seconds")
print(f"⏱️ Minimum total waiting time: {min_wait:.2f} vehicle-seconds")

🔍 Code Breakdown

Inputs:
- lambda_NS, lambda_EW: Define traffic intensities.
- CYCLE_TIME: Fixed to 60 seconds.
Objective Function:
- total_waiting_time(g): Computes total waiting time based on the green time allocated to the NS direction.
- The formula is derived from queueing theory.
Optimization:
- minimize_scalar() searches for the optimal g that minimizes total waiting time within the bounds of 1 to 59 seconds.

📊 Visualizing the Optimization

Let’s visualize how the waiting time changes with different green durations.

green_times = np.linspace(1, CYCLE_TIME - 1, 100)
waiting_times = [total_waiting_time(g) for g in green_times]

plt.figure(figsize=(10, 6))
plt.plot(green_times, waiting_times, label='Total Waiting Time', color='blue')
plt.axvline(optimal_g, color='red', linestyle='--', label=f'Optimal g = {optimal_g:.2f}s')
plt.title('🚦 Traffic Signal Optimization: Total Waiting Time vs Green Time')
plt.xlabel('Green Time for NS (seconds)')
plt.ylabel('Total Waiting Time (vehicle-seconds)')
plt.legend()
plt.grid(True)
plt.show()

📈 Result Analysis

From the graph:

The curve shows a convex shape, confirming the presence of a global minimum.
The red dashed line shows the optimal green time that minimizes waiting.
With asymmetric arrival rates, the optimal green time favors the busier direction (NS in this case).

✅ Conclusion

We successfully:

Modeled a realistic 4-way intersection.
Formulated the total waiting time mathematically.
Used scipy.optimize to find the best green light duration.
Visualized how waiting time depends on signal timing.

This is a great starting point for more advanced systems like:

Multiple intersections (network optimization),
Real-time adaptive traffic signals using reinforcement learning,
Integration with traffic sensors or camera data.

Optimizing Basketball Team Strategy

June 3, 2025

A Data-Driven Approach to Player Rotation

As sports analytics continues to revolutionize how teams make strategic decisions, optimization techniques have become essential tools for coaches and analysts. Today, we’ll explore a practical example of sports strategy optimization using Python, focusing on a basketball team’s player rotation problem.

The Problem: Optimal Player Rotation Strategy

Imagine you’re the analytics consultant for a professional basketball team. The coach wants to optimize player rotations to maximize the team’s offensive efficiency while managing player fatigue throughout a 48-minute game. Each player has different offensive ratings, defensive capabilities, and fatigue rates.

Objective: Maximize total team performance while ensuring no player exceeds their maximum playing time and maintaining at least 5 players on court at all times.

Let’s define our mathematical model:

Decision Variables:

$x_{i,t}$ = 1 if player $i$ plays in time period $t$, 0 otherwise

Objective Function:
$$\max \sum_{i=1}^{n} \sum_{t=1}^{T} (O_i - F_i \cdot t) \cdot x_{i,t}$$

Where:

$O_i$ = offensive rating of player $i$
$F_i$ = fatigue factor of player $i$
$n$ = number of players
$T$ = number of time periods

Constraints:
$$\sum_{t=1}^{T} x_{i,t} \leq M_i \quad \forall i \text{ (max minutes constraint)}$$
$$\sum_{i=1}^{n} x_{i,t} = 5 \quad \forall t \text{ (exactly 5 players on court)}$$

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

# Set up the problem parameters
class BasketballOptimizer:
    def __init__(self):
        # Player data: [Name, Offensive Rating, Defensive Rating, Fatigue Factor, Max Minutes]
        self.players_data = [
            ['LeBron James', 118, 110, 0.8, 42],
            ['Stephen Curry', 125, 105, 1.2, 40],
            ['Kawhi Leonard', 115, 115, 1.0, 35],
            ['Giannis Antetokounmpo', 120, 112, 0.9, 38],
            ['Kevin Durant', 122, 108, 1.1, 36],
            ['Luka Doncic', 116, 104, 1.3, 37],
            ['Jayson Tatum', 114, 109, 1.0, 39],
            ['Joel Embiid', 119, 113, 1.4, 32],
            ['Nikola Jokic', 121, 106, 0.7, 36],
            ['Jimmy Butler', 113, 111, 0.9, 35]
        ]
        
        # Create DataFrame
        self.df = pd.DataFrame(self.players_data, 
                              columns=['Player', 'Off_Rating', 'Def_Rating', 'Fatigue_Factor', 'Max_Minutes'])
        
        # Game parameters
        self.game_duration = 48  # minutes
        self.time_periods = 12   # 4-minute periods
        self.players_on_court = 5
        
    def calculate_performance_matrix(self):
        """Calculate performance for each player in each time period"""
        n_players = len(self.df)
        performance_matrix = np.zeros((n_players, self.time_periods))
        
        for i in range(n_players):
            base_performance = (self.df.iloc[i]['Off_Rating'] + self.df.iloc[i]['Def_Rating']) / 2
            fatigue_factor = self.df.iloc[i]['Fatigue_Factor']
            
            for t in range(self.time_periods):
                # Performance decreases with fatigue over time
                fatigue_penalty = fatigue_factor * (t + 1) * 0.5
                performance_matrix[i, t] = max(base_performance - fatigue_penalty, 50)
                
        return performance_matrix
    
    def solve_optimization(self):
        """Solve the player rotation optimization problem"""
        n_players = len(self.df)
        n_periods = self.time_periods
        
        # Get performance matrix
        performance_matrix = self.calculate_performance_matrix()
        
        # Create decision variables: x[i,t] for player i in period t
        # We'll use a simple greedy approach for demonstration
        rotation_schedule = np.zeros((n_players, n_periods))
        player_minutes = np.zeros(n_players)
        
        # For each time period, select top 5 available players
        for t in range(n_periods):
            # Calculate current performance considering fatigue and availability
            current_scores = []
            for i in range(n_players):
                if player_minutes[i] < self.df.iloc[i]['Max_Minutes']:
                    current_scores.append((performance_matrix[i, t], i))
                else:
                    current_scores.append((0, i))
            
            # Sort by performance and select top 5
            current_scores.sort(reverse=True)
            selected_players = [idx for _, idx in current_scores[:self.players_on_court]]
            
            # Update rotation schedule and minutes
            for player_idx in selected_players:
                rotation_schedule[player_idx, t] = 1
                player_minutes[player_idx] += 4  # 4 minutes per period
        
        return rotation_schedule, performance_matrix
    
    def analyze_results(self, rotation_schedule, performance_matrix):
        """Analyze optimization results"""
        results = {}
        
        # Calculate total minutes per player
        total_minutes = np.sum(rotation_schedule, axis=1) * 4
        
        # Calculate total performance
        total_performance = np.sum(rotation_schedule * performance_matrix)
        
        # Create results dataframe
        results_df = self.df.copy()
        results_df['Minutes_Played'] = total_minutes
        results_df['Utilization'] = (total_minutes / results_df['Max_Minutes'] * 100).round(1)
        
        # Calculate performance contribution
        performance_contribution = np.sum(rotation_schedule * performance_matrix, axis=1)
        results_df['Performance_Contribution'] = performance_contribution.round(1)
        
        results['summary'] = results_df
        results['total_performance'] = total_performance
        results['rotation_matrix'] = rotation_schedule
        
        return results
    
    def create_visualizations(self, results):
        """Create comprehensive visualizations"""
        fig, axes = plt.subplots(2, 2, figsize=(16, 12))
        
        # 1. Player Minutes Distribution
        ax1 = axes[0, 0]
        players = results['summary']['Player']
        minutes = results['summary']['Minutes_Played']
        max_minutes = results['summary']['Max_Minutes']
        
        x_pos = np.arange(len(players))
        bars1 = ax1.bar(x_pos, minutes, alpha=0.7, label='Minutes Played', color='skyblue')
        bars2 = ax1.bar(x_pos, max_minutes, alpha=0.3, label='Max Minutes', color='red')
        
        ax1.set_xlabel('Players')
        ax1.set_ylabel('Minutes')
        ax1.set_title('Player Minutes: Played vs Maximum Available')
        ax1.set_xticks(x_pos)
        ax1.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 2. Performance Contribution
        ax2 = axes[0, 1]
        performance = results['summary']['Performance_Contribution']
        bars = ax2.bar(range(len(players)), performance, color='lightgreen', alpha=0.7)
        
        ax2.set_xlabel('Players')
        ax2.set_ylabel('Performance Contribution')
        ax2.set_title('Individual Player Performance Contribution')
        ax2.set_xticks(range(len(players)))
        ax2.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax2.grid(True, alpha=0.3)
        
        # Add value labels on bars
        for bar, value in zip(bars, performance):
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                    f'{value:.0f}', ha='center', va='bottom')
        
        # 3. Rotation Schedule Heatmap
        ax3 = axes[1, 0]
        rotation_matrix = results['rotation_matrix']
        
        # Create heatmap
        im = ax3.imshow(rotation_matrix, cmap='YlOrRd', aspect='auto')
        ax3.set_xlabel('Time Periods (4-min each)')
        ax3.set_ylabel('Players')
        ax3.set_title('Player Rotation Schedule')
        ax3.set_xticks(range(self.time_periods))
        ax3.set_xticklabels([f'P{i+1}' for i in range(self.time_periods)])
        ax3.set_yticks(range(len(players)))
        ax3.set_yticklabels([name.split()[-1] for name in players])
        
        # Add colorbar
        plt.colorbar(im, ax=ax3, label='Playing (1) / Resting (0)')
        
        # 4. Utilization Rate
        ax4 = axes[1, 1]
        utilization = results['summary']['Utilization']
        colors = ['red' if x > 90 else 'orange' if x > 75 else 'green' for x in utilization]
        
        bars = ax4.bar(range(len(players)), utilization, color=colors, alpha=0.7)
        ax4.axhline(y=100, color='red', linestyle='--', alpha=0.5, label='Max Capacity')
        ax4.axhline(y=75, color='orange', linestyle='--', alpha=0.5, label='High Utilization')
        
        ax4.set_xlabel('Players')
        ax4.set_ylabel('Utilization Rate (%)')
        ax4.set_title('Player Utilization Rate')
        ax4.set_xticks(range(len(players)))
        ax4.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        # Add percentage labels
        for bar, value in zip(bars, utilization):
            ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                    f'{value:.1f}%', ha='center', va='bottom')
        
        plt.tight_layout()
        plt.show()
        
        return fig
    
    def generate_insights(self, results):
        """Generate strategic insights from the optimization"""
        summary = results['summary']
        
        print("=== BASKETBALL TEAM ROTATION OPTIMIZATION RESULTS ===\n")
        
        print("📊 PERFORMANCE SUMMARY:")
        print(f"Total Team Performance Score: {results['total_performance']:.1f}")
        print(f"Average Player Utilization: {summary['Utilization'].mean():.1f}%")
        print(f"Players at >90% Utilization: {sum(summary['Utilization'] > 90)}")
        
        print("\n🏀 TOP PERFORMERS:")
        top_performers = summary.nlargest(3, 'Performance_Contribution')
        for idx, (_, player) in enumerate(top_performers.iterrows()):
            print(f"{idx+1}. {player['Player']}: {player['Performance_Contribution']:.1f} points "
                  f"({player['Minutes_Played']} min, {player['Utilization']:.1f}% utilization)")
        
        print("\n⚠️  STRATEGIC RECOMMENDATIONS:")
        
        # High utilization players
        high_util = summary[summary['Utilization'] > 85]
        if len(high_util) > 0:
            print("• High Utilization Alert:")
            for _, player in high_util.iterrows():
                print(f"  - {player['Player']}: {player['Utilization']:.1f}% utilization - consider rest periods")
        
        # Underutilized players
        low_util = summary[summary['Utilization'] < 60]
        if len(low_util) > 0:
            print("• Underutilized Players:")
            for _, player in low_util.iterrows():
                print(f"  - {player['Player']}: {player['Utilization']:.1f}% utilization - opportunity for more minutes")
        
        # Performance efficiency
        summary['Efficiency'] = summary['Performance_Contribution'] / summary['Minutes_Played']
        top_efficiency = summary.nlargest(3, 'Efficiency')
        print("\n💡 MOST EFFICIENT PLAYERS (Performance per Minute):")
        for idx, (_, player) in enumerate(top_efficiency.iterrows()):
            print(f"{idx+1}. {player['Player']}: {player['Efficiency']:.2f} points/minute")
        
        return summary

# Execute the optimization
def main():
    # Initialize optimizer
    optimizer = BasketballOptimizer()
    
    print("🏀 Basketball Team Rotation Optimization")
    print("=" * 50)
    
    # Display initial team data
    print("\n📋 TEAM ROSTER:")
    print(optimizer.df.to_string(index=False))
    
    # Solve optimization
    print("\n🔄 Solving optimization problem...")
    rotation_schedule, performance_matrix = optimizer.solve_optimization()
    
    # Analyze results
    results = optimizer.analyze_results(rotation_schedule, performance_matrix)
    
    # Generate insights
    insights = optimizer.generate_insights(results)
    
    # Create visualizations
    print("\n📈 Generating visualizations...")
    fig = optimizer.create_visualizations(results)
    
    # Additional analysis: Performance over time
    plt.figure(figsize=(14, 8))
    
    # Plot team performance over time periods
    plt.subplot(2, 1, 1)
    period_performance = np.sum(rotation_schedule * performance_matrix, axis=0)
    periods = [f'Q{(i//3)+1}-{(i%3)*4+4}min' for i in range(optimizer.time_periods)]
    
    plt.plot(periods, period_performance, marker='o', linewidth=2, markersize=8, color='blue')
    plt.title('Team Performance Throughout the Game', fontsize=14, fontweight='bold')
    plt.xlabel('Game Period')
    plt.ylabel('Team Performance Score')
    plt.grid(True, alpha=0.3)
    plt.xticks(rotation=45)
    
    # Plot active players per period
    plt.subplot(2, 1, 2)
    active_players = np.sum(rotation_schedule, axis=0)
    plt.bar(periods, active_players, alpha=0.7, color='lightcoral')
    plt.title('Number of Active Players per Period', fontsize=14, fontweight='bold')
    plt.xlabel('Game Period')
    plt.ylabel('Active Players')
    plt.ylim(0, 6)
    plt.grid(True, alpha=0.3)
    plt.xticks(rotation=45)
    
    # Add horizontal line at 5 players
    plt.axhline(y=5, color='red', linestyle='--', alpha=0.7, label='Required Players')
    plt.legend()
    
    plt.tight_layout()
    plt.show()
    
    return optimizer, results

# Run the main function
if __name__ == "__main__":
    optimizer, results = main()

Code Explanation: Deep Dive into the Basketball Optimization Solution

Let me break down the key components of our basketball rotation optimization system:

1. Data Structure and Initialization

The BasketballOptimizer class begins by defining our player roster with realistic NBA statistics:

Offensive Rating: Points scored per 100 possessions
Defensive Rating: Points allowed per 100 possessions
Fatigue Factor: How quickly performance degrades over time
Max Minutes: Physical limitation for each player

2. Performance Matrix Calculation

The calculate_performance_matrix() method implements our fatigue model:

1 2	fatigue_penalty = fatigue_factor * (t + 1) * 0.5 performance_matrix[i, t] = max(base_performance - fatigue_penalty, 50)

This creates a dynamic performance system where player effectiveness decreases as the game progresses, with different players having different fatigue rates.

3. Optimization Algorithm

Our solve_optimization() method uses a greedy approach that:

Evaluates available players each period based on current performance
Considers minute restrictions and fatigue accumulation
Selects the top 5 performers for each 4-minute period
Updates playing time tracking

4. Results Analysis

The analyze_results() function computes key metrics:

Utilization Rate: $\frac{\text{Minutes Played}}{\text{Max Minutes}} \times 100$
Performance Contribution: $\sum_{t} P_{i,t} \cdot x_{i,t}$
Efficiency: Performance contribution per minute played

🏀 Basketball Team Rotation Optimization
==================================================

📋 TEAM ROSTER:
               Player  Off_Rating  Def_Rating  Fatigue_Factor  Max_Minutes
         LeBron James         118         110             0.8           42
        Stephen Curry         125         105             1.2           40
        Kawhi Leonard         115         115             1.0           35
Giannis Antetokounmpo         120         112             0.9           38
         Kevin Durant         122         108             1.1           36
          Luka Doncic         116         104             1.3           37
         Jayson Tatum         114         109             1.0           39
          Joel Embiid         119         113             1.4           32
         Nikola Jokic         121         106             0.7           36
         Jimmy Butler         113         111             0.9           35

🔄 Solving optimization problem...
=== BASKETBALL TEAM ROTATION OPTIMIZATION RESULTS ===

📊 PERFORMANCE SUMMARY:
Total Team Performance Score: 6698.6
Average Player Utilization: 65.4%
Players at >90% Utilization: 4

🏀 TOP PERFORMERS:
1. Giannis Antetokounmpo: 1135.2 points (40.0 min, 105.3% utilization)
2. Kawhi Leonard: 1012.5 points (36.0 min, 102.9% utilization)
3. Kevin Durant: 1010.2 points (36.0 min, 100.0% utilization)

⚠️  STRATEGIC RECOMMENDATIONS:
• High Utilization Alert:
  - Kawhi Leonard: 102.9% utilization - consider rest periods
  - Giannis Antetokounmpo: 105.3% utilization - consider rest periods
  - Kevin Durant: 100.0% utilization - consider rest periods
  - Joel Embiid: 100.0% utilization - consider rest periods
• Underutilized Players:
  - Luka Doncic: 0.0% utilization - opportunity for more minutes
  - Jayson Tatum: 20.5% utilization - opportunity for more minutes
  - Nikola Jokic: 55.6% utilization - opportunity for more minutes
  - Jimmy Butler: 22.9% utilization - opportunity for more minutes

💡 MOST EFFICIENT PLAYERS (Performance per Minute):
1. Giannis Antetokounmpo: 28.38 points/minute
2. Joel Embiid: 28.17 points/minute
3. Kawhi Leonard: 28.12 points/minute

📈 Generating visualizations...

Visualization Analysis

The comprehensive visualization system provides four critical insights:

Player Minutes Distribution

This chart reveals the balance between player availability and actual usage. Players like LeBron James and Stephen Curry show high utilization rates, indicating they’re core to our strategy.

Performance Contribution

This metric identifies our most impactful players. High performers like Kevin Durant and Giannis Antetokounmpo contribute significantly to overall team success.

Rotation Schedule Heatmap

The heatmap visualization shows our substitution patterns throughout the game. Dark red indicates active periods, revealing strategic rest intervals for key players.

Utilization Rate Analysis

Color-coded bars (green <75%, orange 75-90%, red >90%) help identify potential overuse situations that could lead to fatigue or injury.

Strategic Insights and Recommendations

The optimization reveals several actionable insights:

Load Management: Players with >85% utilization need strategic rest periods
Efficiency Optimization: Some players deliver high performance per minute, suggesting they could handle increased responsibility
Substitution Timing: The performance-over-time graph shows when team effectiveness peaks and dips

Mathematical Foundation

Our optimization problem balances multiple objectives:

Primary Objective:

Subject to constraints:

Resource constraints: $\sum_t x_{i,t} \leq \text{MaxMinutes}_i$
Court requirements: $\sum_i x_{i,t} = 5$
Binary constraints: $x_{i,t} \in {0,1}$

Real-World Applications

This optimization framework extends beyond basketball:

Soccer: Managing player stamina across matches
Hockey: Line changes and shift optimization
Business: Workforce scheduling with skill-based assignments
Military: Resource allocation under constraints

The combination of performance modeling, constraint optimization, and comprehensive visualization provides coaches and analysts with data-driven insights for strategic decision-making. By quantifying the trade-offs between player fatigue, performance, and utilization, teams can maximize their competitive advantage while protecting player health.

This approach demonstrates how mathematical optimization can transform sports strategy from intuition-based decisions to evidence-driven tactical planning.

Detecting Anomalies in Time Series Data

June 2, 2025

A Practical Example with Python

Anomaly detection in time series data is a powerful technique used across industries to identify unusual patterns that deviate from expected behavior. Whether it’s spotting fraudulent transactions, detecting equipment failures, or monitoring system performance, anomaly detection helps us catch issues before they escalate. In this blog post, we’ll walk through a realistic example of detecting anomalies in time series data using Python on Google Colaboratory. We’ll use a synthetic dataset representing hourly server CPU usage, implement a simple yet effective anomaly detection method, and visualize the results with clear explanations.

Problem Statement: Monitoring Server CPU Usage

Imagine you’re a system administrator monitoring the CPU usage of a server over time. The server typically operates within a predictable range, but sudden spikes or drops could indicate issues like overloads, hardware failures, or cyberattacks. Our goal is to detect these anomalies using a statistical approach based on a moving average and standard deviation.

We’ll generate synthetic CPU usage data with some intentional anomalies, apply a z-score-based anomaly detection method, and visualize the results using Matplotlib. Let’s dive into the code and break it down step by step.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

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

# Generate synthetic time series data
dates = [datetime(2025, 6, 1) + timedelta(hours=i) for i in range(168)]  # One week of hourly data
cpu_usage = np.random.normal(loc=50, scale=10, size=len(dates))  # Normal CPU usage ~ N(50, 10)
cpu_usage[50:52] = [90, 95]  # Introduce spike anomalies
cpu_usage[120:122] = [10, 5]  # Introduce drop anomalies

# Create DataFrame
df = pd.DataFrame({'timestamp': dates, 'cpu_usage': cpu_usage})

# Define anomaly detection function
def detect_anomalies(df, window=24, z_threshold=2.5):
    # Calculate rolling mean and standard deviation
    df['rolling_mean'] = df['cpu_usage'].rolling(window=window, center=True).mean()
    df['rolling_std'] = df['cpu_usage'].rolling(window=window, center=True).std()
    
    # Calculate z-score
    df['z_score'] = (df['cpu_usage'] - df['rolling_mean']) / df['rolling_std']
    
    # Identify anomalies
    df['anomaly'] = df['z_score'].abs() > z_threshold
    
    return df

# Apply anomaly detection
df = detect_anomalies(df)

# Plot the results
plt.figure(figsize=(12, 6))
plt.plot(df['timestamp'], df['cpu_usage'], label='CPU Usage', color='blue')
plt.plot(df['timestamp'], df['rolling_mean'], label='Rolling Mean', color='orange', linestyle='--')
plt.fill_between(df['timestamp'], 
                 df['rolling_mean'] - 2.5 * df['rolling_std'], 
                 df['rolling_mean'] + 2.5 * df['rolling_std'], 
                 color='gray', alpha=0.2, label='Normal Range (±2.5σ)')
plt.scatter(df[df['anomaly']]['timestamp'], df[df['anomaly']]['cpu_usage'], 
            color='red', label='Anomalies', s=100, marker='o')
plt.title('Server CPU Usage with Anomaly Detection')
plt.xlabel('Timestamp')
plt.ylabel('CPU Usage (%)')
plt.legend()
plt.grid(True)
plt.xticks(rotation=45)
plt.tight_layout()

# Save the plot
plt.savefig('cpu_usage_anomaly_detection.png')

Code Explanation

Let’s break down the Python code to understand how it detects anomalies and visualizes the results.

Import Libraries:
- numpy and pandas are used for numerical operations and data handling.
- matplotlib.pyplot is used for plotting.
- datetime and timedelta help create a time series of hourly timestamps.
Generate Synthetic Data:
- We create a week’s worth of hourly data (168 hours) starting from June 1, 2025.
- CPU usage is modeled as a normal distribution with mean $(\mu = 50)$ and standard deviation $(\sigma = 10)$, i.e., $( \text{cpu_usage} \sim \mathcal{N}(50, 10) )$.
- We introduce deliberate anomalies:
  - At hours 50–51, we set CPU usage to 90% and 95% to simulate spikes.
  - At hours 120–121, we set CPU usage to 10% and 5% to simulate drops.
- The data is stored in a Pandas DataFrame with timestamp and cpu_usage columns.
Anomaly Detection Function:
- The detect_anomalies function uses a rolling window approach to compute the moving average and standard deviation.
- Rolling Mean: We calculate the mean over a window of 24 hours (centered) to capture the typical CPU usage trend:
- Rolling Standard Deviation: We compute the standard deviation over the same window to measure variability:

Z-Score: For each data point, we calculate the z-score to measure how far it deviates from the mean in terms of standard deviations:
Anomaly Detection: A point is flagged as an anomaly if , where we set This corresponds to points outside approximately 99% of the data under a normal distribution.

Plotting:
- We create a plot with:
  - The CPU usage time series (blue line).
  - The rolling mean (orange dashed line).
  - A shaded band representing the “normal” range $(( \text{rolling_mean} \pm 2.5 \cdot \text{rolling_std} )$).
  - Red dots marking detected anomalies.
- The plot is saved as cpu_usage_anomaly_detection.png.

Visualizing and Interpreting the Results

The plot generated by the code provides a clear view of the CPU usage over time and highlights anomalies. Here’s what you’ll see:

CPU Usage (Blue Line): This shows the hourly CPU usage, fluctuating around 50% with some visible spikes and drops.
Rolling Mean (Orange Dashed Line): This represents the average CPU usage over a 24-hour window, smoothing out short-term fluctuations.
Normal Range (Gray Band): The shaded area shows the expected range of CPU usage $(( \mu \pm 2.5\sigma )$). Most data points should fall within this band.
Anomalies (Red Dots): Points outside the gray band (where $( |z_t| > 2.5 )$) are marked as anomalies. You’ll notice red dots at hours 50–51 (spikes to 90% and 95%) and hours 120–121 (drops to 10% and 5%).

Interpretation:

The spikes at hours 50–51 could indicate a server overload, perhaps due to a sudden surge in traffic or a resource-intensive process.
The drops at hours 120–121 might suggest a system failure or a period of inactivity, warranting investigation.
The z-score threshold of 2.5 ensures we only flag significant deviations, reducing false positives while catching critical anomalies.

Why This Approach Works

Using a rolling mean and standard deviation is a simple yet effective method for anomaly detection in time series data. It’s robust for data with stable patterns and can adapt to gradual changes in the mean or variance. The z-score provides a standardized way to measure deviations, making it easy to set a threshold (e.g., 2.5σ) based on the desired sensitivity.

However, this method assumes the data is roughly normally distributed within the window. For more complex time series with seasonality or trends, you might consider advanced techniques like ARIMA, Prophet, or machine learning models like Isolation Forest. For our server CPU usage scenario, this approach is practical and computationally efficient.

Running the Code

To try this yourself, copy the code into a Google Colaboratory notebook and run it. The plot will be saved as cpu_usage_anomaly_detection.png in your Colab environment. You can download it to inspect the results or modify parameters like the window size (window) or z-score threshold (z_threshold) to experiment with sensitivity.

This example demonstrates how anomaly detection can be applied to real-world problems with minimal code. Whether you’re monitoring servers, financial transactions, or sensor data, this approach is a great starting point for identifying outliers in time series data.

Happy coding, and stay vigilant for those anomalies! 🚨

Estimating the Position of a Noisy GPS

June 1, 2025

🚗 A Kalman Filter Example in Python

In this post, we’ll explore a real-world application of the Kalman Filter using Python: tracking a car’s position using noisy GPS measurements.
GPS sensors are often subject to random noise, so Kalman filters are widely used in navigation systems to smooth and predict location.

Let’s dive in with a concrete scenario.

🧠 Problem Statement

Imagine a car driving along a straight road at a constant speed.
The GPS gives us the position every second, but with noise.
We want to estimate the true position and velocity of the car over time using the Kalman Filter.

🧮 Mathematical Model

We’ll track two state variables:

$x$: position
$v$: velocity

The state vector is:

$$
\mathbf{x} =
\begin{bmatrix}
x \
v
\end{bmatrix}
$$

We model the state transition as:

Where $\Delta t$ is the time step, and $\mathbf{w}_k$ is process noise.

The measurement (GPS) only gives the position:

$$
\mathbf{z}_k =
\begin{bmatrix}
1 & 0
\end{bmatrix}
\mathbf{x}_k + \mathbf{v}_k
$$

Where $\mathbf{v}_k$ is measurement noise.

🧪 Implementation in Python

Let’s implement this in Python using NumPy and Matplotlib.

import numpy as np
import matplotlib.pyplot as plt

# Time parameters
dt = 1.0  # 1 second between measurements
N = 50    # number of time steps

# True position and velocity
true_velocity = 1.0
true_position = np.arange(N) * true_velocity

# Simulated noisy GPS measurements
np.random.seed(42)
measurement_noise_std = 2.0
measurements = true_position + np.random.normal(0, measurement_noise_std, size=N)

# Kalman Filter variables
x = np.array([[0], [0]])  # Initial state estimate: position=0, velocity=0
P = np.eye(2) * 500        # Initial uncertainty

F = np.array([[1, dt],
              [0, 1]])     # State transition matrix
H = np.array([[1, 0]])     # Measurement matrix

R = np.array([[measurement_noise_std**2]])  # Measurement noise covariance
Q = np.array([[0.1, 0.0],                   # Process noise covariance
              [0.0, 0.1]])

# Storage for estimates
estimated_positions = []
estimated_velocities = []

for z in measurements:
    # Prediction
    x = F @ x
    P = F @ P @ F.T + Q

    # Measurement update
    y = np.array([[z]]) - (H @ x)
    S = H @ P @ H.T + R
    K = P @ H.T @ np.linalg.inv(S)
    x = x + K @ y
    P = (np.eye(2) - K @ H) @ P

    # Store results
    estimated_positions.append(x[0, 0])
    estimated_velocities.append(x[1, 0])

🔍 Code Explanation

F: State transition matrix incorporates motion ($x_{k+1} = x_k + v_k \Delta t$).
H: Measurement matrix — we only observe position.
Q: Models uncertainty in the system’s dynamics (e.g., acceleration not modeled).
R: Models measurement noise from the GPS.

Each iteration involves:

Prediction: Estimate next state based on current velocity.
Update: Correct estimate using noisy GPS measurement.

📊 Visualization of Results

Let’s plot the noisy measurements, true position, and Kalman filter estimates:

plt.figure(figsize=(12, 6))
plt.plot(true_position, label="True Position", linewidth=2)
plt.plot(measurements, label="Noisy GPS Measurements", linestyle='dotted')
plt.plot(estimated_positions, label="Kalman Filter Estimate", linewidth=2)
plt.xlabel("Time step")
plt.ylabel("Position")
plt.legend()
plt.title("Kalman Filter: GPS Position Tracking")
plt.grid(True)
plt.show()

🧭 Velocity Estimation

We can also plot how the velocity estimate converges:

plt.figure(figsize=(12, 4))
plt.plot(estimated_velocities, label="Estimated Velocity", color="orange")
plt.axhline(true_velocity, color='green', linestyle='--', label="True Velocity")
plt.xlabel("Time step")
plt.ylabel("Velocity")
plt.legend()
plt.title("Kalman Filter: Velocity Estimation")
plt.grid(True)
plt.show()

✅ Interpretation of Results

Initially, the Kalman filter struggles due to high initial uncertainty.
Over time, it learns both the true position and velocity.
Despite noisy GPS, the estimates become stable and accurate.
This method is powerful for real-time applications like self-driving cars, drones, and robotics.

🔚 Conclusion

The Kalman filter is an elegant tool to combine predictions and noisy measurements.
With just a few lines of linear algebra, we’ve created a robust estimator for tracking motion in the presence of noise.

Principal Component Analysis

May 31, 2025

A Practical Guide with Python

Today, we’ll dive into Principal Component Analysis (PCA) using a real-world dataset - the famous Wine dataset. This example will demonstrate how PCA can help us understand high-dimensional data by reducing its complexity while preserving the most important information.

The Problem

Imagine you’re a wine researcher with measurements of 13 different chemical properties for 178 wine samples from three different cultivars. With 13 dimensions, it’s impossible to visualize the data effectively. PCA will help us reduce these dimensions while maintaining the essential characteristics that distinguish different wine types.

Mathematical Foundation

PCA finds the directions (principal components) along which the data varies the most. Mathematically, we:

Standardize the data: $X_{standardized} = \frac{X - \mu}{\sigma}$
Compute the covariance matrix: $C = \frac{1}{n-1}X^TX$
Find eigenvalues and eigenvectors: $C\mathbf{v} = \lambda\mathbf{v}$
Transform the data: $Y = X \cdot \mathbf{V}$

Where $\mathbf{V}$ contains the eigenvectors (principal components) and $\lambda$ are the eigenvalues representing the variance explained by each component.

# Import necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D

# Set style for better visualizations
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

# Load the Wine dataset
wine_data = load_wine()
X = wine_data.data
y = wine_data.target
feature_names = wine_data.feature_names
target_names = wine_data.target_names

print("Dataset Information:")
print(f"Number of samples: {X.shape[0]}")
print(f"Number of features: {X.shape[1]}")
print(f"Target classes: {target_names}")
print(f"Class distribution: {np.bincount(y)}")

# Convert to DataFrame for easier handling
df = pd.DataFrame(X, columns=feature_names)
df['target'] = y
df['target_name'] = [target_names[i] for i in y]

print("\nFirst 5 rows of the dataset:")
print(df.head())

# Display basic statistics
print("\nDataset Statistics:")
print(df.describe())

# Step 1: Data Standardization
print("\n" + "="*50)
print("STEP 1: DATA STANDARDIZATION")
print("="*50)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("Original data statistics:")
print(f"Mean: {np.mean(X, axis=0)[:3]}... (showing first 3 features)")
print(f"Std:  {np.std(X, axis=0)[:3]}... (showing first 3 features)")

print("\nStandardized data statistics:")
print(f"Mean: {np.mean(X_scaled, axis=0)[:3]}... (should be ~0)")
print(f"Std:  {np.std(X_scaled, axis=0)[:3]}... (should be ~1)")

# Step 2: Apply PCA
print("\n" + "="*50)
print("STEP 2: PRINCIPAL COMPONENT ANALYSIS")
print("="*50)

# Perform PCA with all components first to analyze explained variance
pca_full = PCA()
X_pca_full = pca_full.fit_transform(X_scaled)

# Calculate cumulative explained variance
explained_variance_ratio = pca_full.explained_variance_ratio_
cumulative_variance = np.cumsum(explained_variance_ratio)

print("Explained Variance by each Principal Component:")
for i, (var, cum_var) in enumerate(zip(explained_variance_ratio, cumulative_variance)):
    if i < 10:  # Show first 10 components
        print(f"PC{i+1}: {var:.4f} ({var*100:.2f}%) - Cumulative: {cum_var:.4f} ({cum_var*100:.2f}%)")

# Determine optimal number of components (e.g., 95% variance)
n_components_95 = np.argmax(cumulative_variance >= 0.95) + 1
n_components_90 = np.argmax(cumulative_variance >= 0.90) + 1

print(f"\nComponents needed for 90% variance: {n_components_90}")
print(f"Components needed for 95% variance: {n_components_95}")

# Apply PCA with reduced dimensions (let's use 3 components for visualization)
pca_3d = PCA(n_components=3)
X_pca_3d = pca_3d.fit_transform(X_scaled)

print(f"\nExplained variance with 3 components: {sum(pca_3d.explained_variance_ratio_):.4f} ({sum(pca_3d.explained_variance_ratio_)*100:.2f}%)")

# Step 3: Analyze Principal Components
print("\n" + "="*50)
print("STEP 3: PRINCIPAL COMPONENTS ANALYSIS")
print("="*50)

# Create DataFrame for PC loadings
pc_loadings = pd.DataFrame(
    pca_3d.components_.T,
    columns=['PC1', 'PC2', 'PC3'],
    index=feature_names
)

print("Principal Component Loadings (Top 5 features for each PC):")
for pc in ['PC1', 'PC2', 'PC3']:
    print(f"\n{pc} - Top contributors:")
    top_features = pc_loadings[pc].abs().sort_values(ascending=False).head(5)
    for feature, loading in top_features.items():
        direction = "+" if pc_loadings.loc[feature, pc] > 0 else "-"
        print(f"  {direction} {feature}: {abs(loading):.3f}")

# Step 4: Visualization
print("\n" + "="*50)
print("STEP 4: CREATING VISUALIZATIONS")
print("="*50)

# Create comprehensive visualization
fig = plt.figure(figsize=(20, 15))

# 1. Explained Variance Plot
ax1 = plt.subplot(2, 3, 1)
plt.bar(range(1, 14), explained_variance_ratio, alpha=0.7, label='Individual')
plt.plot(range(1, 14), cumulative_variance, 'ro-', label='Cumulative')
plt.xlabel('Principal Component')
plt.ylabel('Explained Variance Ratio')
plt.title('Explained Variance by Principal Components')
plt.legend()
plt.grid(True, alpha=0.3)

# Add annotation for 95% threshold
plt.axhline(y=0.95, color='red', linestyle='--', alpha=0.7)
plt.text(8, 0.96, '95% threshold', fontsize=10)

# 2. 2D PCA Scatter Plot (PC1 vs PC2)
ax2 = plt.subplot(2, 3, 2)
colors = ['red', 'blue', 'green']
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    plt.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 1], 
               c=color, label=target_class, alpha=0.7, s=50)

plt.xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%} variance)')
plt.ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%} variance)')
plt.title('PCA: First Two Principal Components')
plt.legend()
plt.grid(True, alpha=0.3)

# 3. 3D PCA Scatter Plot
ax3 = plt.subplot(2, 3, 3, projection='3d')
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    ax3.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 1], X_pca_3d[mask, 2],
               c=color, label=target_class, alpha=0.7, s=30)

ax3.set_xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%})')
ax3.set_ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%})')
ax3.set_zlabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%})')
ax3.set_title('3D PCA Visualization')
ax3.legend()

# 4. Feature Contributions Heatmap
ax4 = plt.subplot(2, 3, 4)
sns.heatmap(pc_loadings.T, annot=False, cmap='RdBu_r', center=0,
            xticklabels=[name.replace('_', ' ').title() for name in feature_names],
            yticklabels=['PC1', 'PC2', 'PC3'])
plt.title('Principal Component Loadings Heatmap')
plt.xlabel('Features')
plt.ylabel('Principal Components')
plt.xticks(rotation=45, ha='right')

# 5. PC1 vs PC3
ax5 = plt.subplot(2, 3, 5)
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    plt.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 2], 
               c=color, label=target_class, alpha=0.7, s=50)

plt.xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%} variance)')
plt.ylabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%} variance)')
plt.title('PCA: PC1 vs PC3')
plt.legend()
plt.grid(True, alpha=0.3)

# 6. PC2 vs PC3
ax6 = plt.subplot(2, 3, 6)
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    plt.scatter(X_pca_3d[mask, 1], X_pca_3d[mask, 2], 
               c=color, label=target_class, alpha=0.7, s=50)

plt.xlabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%} variance)')
plt.ylabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%} variance)')
plt.title('PCA: PC2 vs PC3')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional Analysis: Biplot
print("\nCreating Biplot for detailed feature analysis...")
fig, ax = plt.subplots(figsize=(12, 10))

# Plot data points
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    ax.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 1], 
              c=color, label=target_class, alpha=0.6, s=50)

# Plot feature vectors
feature_vectors = pca_3d.components_[:2].T * 3  # Scale for visibility
for i, (feature, vector) in enumerate(zip(feature_names, feature_vectors)):
    ax.arrow(0, 0, vector[0], vector[1], head_width=0.1, head_length=0.1, 
             fc='black', ec='black', alpha=0.7)
    ax.text(vector[0]*1.1, vector[1]*1.1, feature.replace('_', ' ').title(), 
            fontsize=8, ha='center', va='center')

ax.set_xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%} variance)')
ax.set_ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%} variance)')
ax.set_title('PCA Biplot: Data Points and Feature Contributions')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)

plt.tight_layout()
plt.show()

# Step 5: Reconstruction Analysis
print("\n" + "="*50)
print("STEP 5: RECONSTRUCTION ANALYSIS")
print("="*50)

# Reconstruct data using different numbers of components
n_components_list = [1, 2, 3, 5, 7, 10, 13]
reconstruction_errors = []

for n_comp in n_components_list:
    pca_temp = PCA(n_components=n_comp)
    X_transformed = pca_temp.fit_transform(X_scaled)
    X_reconstructed = pca_temp.inverse_transform(X_transformed)
    
    # Calculate reconstruction error (MSE)
    mse = np.mean((X_scaled - X_reconstructed) ** 2)
    reconstruction_errors.append(mse)
    
    if n_comp <= 5:
        print(f"Components: {n_comp:2d}, Explained Variance: {sum(pca_temp.explained_variance_ratio_):.4f}, Reconstruction Error: {mse:.6f}")

# Plot reconstruction error
plt.figure(figsize=(10, 6))
plt.plot(n_components_list, reconstruction_errors, 'bo-', linewidth=2, markersize=8)
plt.xlabel('Number of Principal Components')
plt.ylabel('Reconstruction Error (MSE)')
plt.title('Reconstruction Error vs Number of Components')
plt.grid(True, alpha=0.3)
plt.yscale('log')

# Add annotations
for i, (n_comp, error) in enumerate(zip(n_components_list, reconstruction_errors)):
    if n_comp in [2, 3, 5]:
        plt.annotate(f'{error:.4f}', (n_comp, error), 
                    textcoords="offset points", xytext=(0,10), ha='center')

plt.tight_layout()
plt.show()

# Final Summary
print("\n" + "="*50)
print("ANALYSIS SUMMARY")
print("="*50)
print(f"• Original dataset: {X.shape[0]} samples × {X.shape[1]} features")
print(f"• First 3 PCs explain {sum(pca_3d.explained_variance_ratio_)*100:.1f}% of total variance")
print(f"• PC1 captures {pca_3d.explained_variance_ratio_[0]*100:.1f}% (likely represents overall wine quality/intensity)")
print(f"• PC2 captures {pca_3d.explained_variance_ratio_[1]*100:.1f}% (likely represents color-related properties)")
print(f"• PC3 captures {pca_3d.explained_variance_ratio_[2]*100:.1f}% (likely represents acidity/pH balance)")
print(f"• The three wine classes show good separation in the PC space")
print(f"• Dimensionality reduction from 13D to 3D maintains {sum(pca_3d.explained_variance_ratio_)*100:.1f}% of information")

Code Explanation

Let me break down the key components of this PCA implementation:

1. Data Loading and Preparation

The code starts by loading the Wine dataset from scikit-learn, which contains 13 chemical measurements for 178 wine samples from three different cultivars. We examine the basic statistics to understand our data structure.

2. Data Standardization

1 2	scaler = StandardScaler() X_scaled = scaler.fit_transform(X)

This step is crucial because PCA is sensitive to the scale of features. The formula $X_{standardized} = \frac{X - \mu}{\sigma}$ ensures all features have mean = 0 and standard deviation = 1.

3. PCA Implementation

The code implements PCA in two phases:

Full PCA: Analyzes all 13 components to understand the variance distribution
Reduced PCA: Focuses on the first 3 components for visualization

The eigenvalue decomposition finds the directions of maximum variance: $C\mathbf{v} = \lambda\mathbf{v}$

4. Component Analysis

The loadings matrix shows how much each original feature contributes to each principal component. High absolute values indicate strong influence on that component’s direction.

5. Comprehensive Visualization

The visualization includes six different plots:

Explained Variance: Shows how much information each PC captures
2D/3D Scatter Plots: Display the transformed data in reduced dimensions
Heatmap: Visualizes feature contributions to each PC
Biplot: Combines data points with feature vectors for interpretation

6. Reconstruction Analysis

This measures how well we can recreate the original data using fewer dimensions. The reconstruction error formula is:
$$\text{MSE} = \frac{1}{n \times p} \sum_{i=1}^{n} \sum_{j=1}^{p} (X_{ij} - \hat{X}_{ij})^2$$

Results

Dataset Information:
Number of samples: 178
Number of features: 13
Target classes: ['class_0' 'class_1' 'class_2']
Class distribution: [59 71 48]

First 5 rows of the dataset:
   alcohol  malic_acid   ash  alcalinity_of_ash  magnesium  total_phenols  \
0    14.23        1.71  2.43               15.6      127.0           2.80   
1    13.20        1.78  2.14               11.2      100.0           2.65   
2    13.16        2.36  2.67               18.6      101.0           2.80   
3    14.37        1.95  2.50               16.8      113.0           3.85   
4    13.24        2.59  2.87               21.0      118.0           2.80   

   flavanoids  nonflavanoid_phenols  proanthocyanins  color_intensity   hue  \
0        3.06                  0.28             2.29             5.64  1.04   
1        2.76                  0.26             1.28             4.38  1.05   
2        3.24                  0.30             2.81             5.68  1.03   
3        3.49                  0.24             2.18             7.80  0.86   
4        2.69                  0.39             1.82             4.32  1.04   

   od280/od315_of_diluted_wines  proline  target target_name  
0                          3.92   1065.0       0     class_0  
1                          3.40   1050.0       0     class_0  
2                          3.17   1185.0       0     class_0  
3                          3.45   1480.0       0     class_0  
4                          2.93    735.0       0     class_0  

Dataset Statistics:
          alcohol  malic_acid         ash  alcalinity_of_ash   magnesium  \
count  178.000000  178.000000  178.000000         178.000000  178.000000   
mean    13.000618    2.336348    2.366517          19.494944   99.741573   
std      0.811827    1.117146    0.274344           3.339564   14.282484   
min     11.030000    0.740000    1.360000          10.600000   70.000000   
25%     12.362500    1.602500    2.210000          17.200000   88.000000   
50%     13.050000    1.865000    2.360000          19.500000   98.000000   
75%     13.677500    3.082500    2.557500          21.500000  107.000000   
max     14.830000    5.800000    3.230000          30.000000  162.000000   

       total_phenols  flavanoids  nonflavanoid_phenols  proanthocyanins  \
count     178.000000  178.000000            178.000000       178.000000   
mean        2.295112    2.029270              0.361854         1.590899   
std         0.625851    0.998859              0.124453         0.572359   
min         0.980000    0.340000              0.130000         0.410000   
25%         1.742500    1.205000              0.270000         1.250000   
50%         2.355000    2.135000              0.340000         1.555000   
75%         2.800000    2.875000              0.437500         1.950000   
max         3.880000    5.080000              0.660000         3.580000   

       color_intensity         hue  od280/od315_of_diluted_wines      proline  \
count       178.000000  178.000000                    178.000000   178.000000   
mean          5.058090    0.957449                      2.611685   746.893258   
std           2.318286    0.228572                      0.709990   314.907474   
min           1.280000    0.480000                      1.270000   278.000000   
25%           3.220000    0.782500                      1.937500   500.500000   
50%           4.690000    0.965000                      2.780000   673.500000   
75%           6.200000    1.120000                      3.170000   985.000000   
max          13.000000    1.710000                      4.000000  1680.000000   

           target  
count  178.000000  
mean     0.938202  
std      0.775035  
min      0.000000  
25%      0.000000  
50%      1.000000  
75%      2.000000  
max      2.000000  

==================================================
STEP 1: DATA STANDARDIZATION
==================================================
Original data statistics:
Mean: [13.00061798  2.33634831  2.36651685]... (showing first 3 features)
Std:  [0.80954291 1.11400363 0.27357229]... (showing first 3 features)

Standardized data statistics:
Mean: [ 7.84141790e-15  2.44498554e-16 -4.05917497e-15]... (should be ~0)
Std:  [1. 1. 1.]... (should be ~1)

==================================================
STEP 2: PRINCIPAL COMPONENT ANALYSIS
==================================================
Explained Variance by each Principal Component:
PC1: 0.3620 (36.20%) - Cumulative: 0.3620 (36.20%)
PC2: 0.1921 (19.21%) - Cumulative: 0.5541 (55.41%)
PC3: 0.1112 (11.12%) - Cumulative: 0.6653 (66.53%)
PC4: 0.0707 (7.07%) - Cumulative: 0.7360 (73.60%)
PC5: 0.0656 (6.56%) - Cumulative: 0.8016 (80.16%)
PC6: 0.0494 (4.94%) - Cumulative: 0.8510 (85.10%)
PC7: 0.0424 (4.24%) - Cumulative: 0.8934 (89.34%)
PC8: 0.0268 (2.68%) - Cumulative: 0.9202 (92.02%)
PC9: 0.0222 (2.22%) - Cumulative: 0.9424 (94.24%)
PC10: 0.0193 (1.93%) - Cumulative: 0.9617 (96.17%)

Components needed for 90% variance: 8
Components needed for 95% variance: 10

Explained variance with 3 components: 0.6653 (66.53%)

==================================================
STEP 3: PRINCIPAL COMPONENTS ANALYSIS
==================================================
Principal Component Loadings (Top 5 features for each PC):

PC1 - Top contributors:
  + flavanoids: 0.423
  + total_phenols: 0.395
  + od280/od315_of_diluted_wines: 0.376
  + proanthocyanins: 0.313
  - nonflavanoid_phenols: 0.299

PC2 - Top contributors:
  + color_intensity: 0.530
  + alcohol: 0.484
  + proline: 0.365
  + ash: 0.316
  + magnesium: 0.300

PC3 - Top contributors:
  + ash: 0.626
  + alcalinity_of_ash: 0.612
  - alcohol: 0.207
  + nonflavanoid_phenols: 0.170
  + od280/od315_of_diluted_wines: 0.166

==================================================
STEP 4: CREATING VISUALIZATIONS
==================================================

Creating Biplot for detailed feature analysis...

==================================================
STEP 5: RECONSTRUCTION ANALYSIS
==================================================
Components:  1, Explained Variance: 0.3620, Reconstruction Error: 0.638012
Components:  2, Explained Variance: 0.5541, Reconstruction Error: 0.445937
Components:  3, Explained Variance: 0.6653, Reconstruction Error: 0.334700
Components:  5, Explained Variance: 0.8016, Reconstruction Error: 0.198377

==================================================
ANALYSIS SUMMARY
==================================================
• Original dataset: 178 samples × 13 features
• First 3 PCs explain 66.5% of total variance
• PC1 captures 36.2% (likely represents overall wine quality/intensity)
• PC2 captures 19.2% (likely represents color-related properties)
• PC3 captures 11.1% (likely represents acidity/pH balance)
• The three wine classes show good separation in the PC space
• Dimensionality reduction from 13D to 3D maintains 66.5% of information

Key Insights from the Results

Dimensionality Reduction: The first 3 principal components typically capture around 66-70% of the total variance, allowing us to visualize 13-dimensional data effectively.
Class Separation: The three wine cultivars show distinct clustering in the PCA space, indicating that the chemical properties effectively distinguish between wine types.
Feature Interpretation:
- PC1 often relates to overall wine intensity/quality
- PC2 typically captures color-related properties
- PC3 usually represents acidity and pH balance
Practical Application: This analysis helps wine researchers understand which chemical properties are most important for distinguishing wine types and can guide quality control processes.

The PCA transformation formula $Y = X \cdot \mathbf{V}$ projects our 13-dimensional wine data into a 3-dimensional space while preserving the most significant patterns, making complex data interpretable and visualizable.

Understanding Ridge Regression with a Practical Example in Python

May 30, 2025

Today, we’re diving into Ridge Regression, a powerful technique for handling regression problems, especially when dealing with multicollinearity or overfitting.
We’ll walk through a concrete example, implement it in Python, and visualize the results to make sense of it all.

Let’s get started!

What is Ridge Regression?

Ridge Regression is a type of linear regression that includes a regularization term to prevent overfitting.
The standard linear regression model minimizes the sum of squared residuals:

$$
\text{Minimize} \sum_{i=1}^n (y_i - \hat{y}_i)^2
$$

where $( y_i )$ is the actual value, and is the predicted value.

Ridge Regression adds an L2 penalty term to the cost function:

Here, $( \lambda )$ (the regularization parameter) controls the strength of the penalty.
A larger $( \lambda )$ shrinks the coefficients $( \beta_j )$ toward zero, reducing model complexity and preventing overfitting.

The Example Problem

Let’s create a synthetic dataset to demonstrate Ridge Regression.
Imagine we’re predicting house prices ($( y )$) based on three features: house size (in square feet), number of bedrooms, and age of the house (in years).
These features might be correlated (e.g., larger houses often have more bedrooms), making Ridge Regression a great choice to stabilize the model.

We’ll generate a dataset with 100 samples, add some noise, and include multicollinearity between features.
Then, we’ll fit both a standard Linear Regression model and a Ridge Regression model to compare their performance.

Python Implementation

Below is the complete Python code to generate the data, fit the models, and visualize the results.
You can run this directly in Google Colaboratory.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler

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

# Generate synthetic dataset
n_samples = 100
house_size = np.random.normal(2000, 500, n_samples)  # Size in sq ft
bedrooms = house_size / 500 + np.random.normal(0, 0.5, n_samples)  # Correlated with size
age = np.random.normal(20, 10, n_samples)  # Age in years
X = np.column_stack((house_size, bedrooms, age))

# True coefficients: price = 100 * size + 50 * bedrooms - 20 * age + noise
true_coefficients = [100, 50, -20]
y = X @ true_coefficients + np.random.normal(0, 10000, n_samples)

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Fit Linear Regression
lr = LinearRegression()
lr.fit(X_scaled, y)

# Fit Ridge Regression with different alpha values
alphas = [0.1, 1.0, 10.0]
ridge_models = {alpha: Ridge(alpha=alpha).fit(X_scaled, y) for alpha in alphas}

# Predictions
y_pred_lr = lr.predict(X_scaled)
y_pred_ridge = {alpha: model.predict(X_scaled) for alpha, model in ridge_models.items()}

# Calculate MSE
mse_lr = mean_squared_error(y, y_pred_lr)
mse_ridge = {alpha: mean_squared_error(y, y_pred_ridge[alpha]) for alpha in alphas}

# Plotting
plt.figure(figsize=(12, 6))

# Plot 1: Actual vs Predicted for Linear Regression
plt.subplot(1, 2, 1)
plt.scatter(y, y_pred_lr, color='blue', alpha=0.5, label='Predicted')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2, label='Ideal')
plt.xlabel('Actual Price ($)')
plt.ylabel('Predicted Price ($)')
plt.title(f'Linear Regression\nMSE: {mse_lr:.2f}')
plt.legend()

# Plot 2: Actual vs Predicted for Ridge Regression (alpha=1.0)
plt.subplot(1, 2, 2)
plt.scatter(y, y_pred_ridge[1.0], color='green', alpha=0.5, label='Predicted')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2, label='Ideal')
plt.xlabel('Actual Price ($)')
plt.ylabel('Predicted Price ($)')
plt.title(f'Ridge Regression (α=1.0)\nMSE: {mse_ridge[1.0]:.2f}')
plt.legend()

plt.tight_layout()
plt.savefig('ridge_regression_plot.png')
plt.show()

# Print coefficients
print("Linear Regression Coefficients:", lr.coef_)
print("Ridge Regression Coefficients:")
for alpha, model in ridge_models.items():
    print(f"Alpha = {alpha}: {model.coef_}")

Code Explanation

Let’s break down the code step by step:

Import Libraries:
- We use numpy for numerical operations, pandas (though minimally here), and matplotlib.pyplot for plotting.
- From sklearn, we import LinearRegression for the baseline model, Ridge for Ridge Regression, mean_squared_error to evaluate performance, and StandardScaler to standardize features.
Generate Synthetic Data:
- We set a random seed (np.random.seed(42)) for reproducibility.
- We create 100 samples:
  - house_size: Normally distributed around 2000 sq ft with a standard deviation of 500.
  - bedrooms: Correlated with house_size (divided by 500 with added noise) to simulate multicollinearity.
  - age: Normally distributed around 20 years with a standard deviation of 10.
- Features are stacked into a matrix X using np.column_stack.
- The target y (house price) is computed as a linear combination of features with true coefficients [100, 50, -20] plus Gaussian noise.
Feature Scaling:
- We use StandardScaler to standardize X (mean = 0, variance = 1). This is crucial for Ridge Regression because the L2 penalty is sensitive to feature scales.
Model Fitting:
- A LinearRegression model is fitted to the scaled data.
- Ridge Regression models are fitted with three different $( \lambda )$ values (alpha = 0.1, 1.0, 10.0) to explore the effect of regularization strength.
Predictions and Evaluation:
- We generate predictions for both models.
- Mean Squared Error (MSE) is calculated to quantify prediction accuracy.
Visualization:
- Two scatter plots are created:
  - Left: Actual vs. predicted prices for Linear Regression.
  - Right: Actual vs. predicted prices for Ridge Regression with $( \alpha = 1.0 )$.
- A red dashed line represents the ideal case (perfect predictions).
- MSE is displayed in the plot titles.
- The plot is saved as ridge_regression_plot.png and displayed.
Print Coefficients:
- We print the coefficients of both models to compare how Ridge Regression shrinks them compared to Linear Regression.

Visualizing and Interpreting the Results

Running the code produces a plot with two subplots:

Left Subplot (Linear Regression):
- Shows actual house prices (x-axis) vs. predicted prices (y-axis).
- Points close to the red dashed line indicate accurate predictions.
- The MSE (e.g., ~100,000,000) reflects the model’s error on the training data.
Right Subplot (Ridge Regression, $( \alpha = 1.0 )$):
- Similarly shows actual vs. predicted prices.
- The MSE is typically close to or slightly better than Linear Regression, indicating that regularization helps stabilize predictions.

The printed coefficients reveal the effect of regularization:

Linear Regression Coefficients: These may be large due to multicollinearity between house_size and bedrooms.
Ridge Regression Coefficients: As $( \alpha )$ increases, coefficients shrink toward zero, reducing the model’s sensitivity to correlated features.

For example, you might see output like:

Linear Regression Coefficients: [44269.64703873  -687.42515039    74.44582334]
Ridge Regression Coefficients:
Alpha = 0.1: [44081.22365278  -524.05649958    83.97555061]
Alpha = 1.0: [42501.98691588   831.74436073   165.15154521]
Alpha = 10.0: [33190.91186375  8068.96266086   706.04219055]

Notice how Ridge Regression reduces the magnitude of coefficients as $( \alpha )$ increases, mitigating overfitting.

Why Ridge Regression Shines

In our example, the correlation between house_size and bedrooms simulates a real-world scenario where features are not independent.
Linear Regression might overfit by assigning inflated coefficients to correlated features.
Ridge Regression, by contrast, balances the trade-off between fitting the data and keeping coefficients small, leading to a more robust model.

Conclusion

Ridge Regression is a fantastic tool for regression tasks with correlated features or potential overfitting.
Our example showed how to implement it in Python, compare it to Linear Regression, and visualize the results.
The scatter plots and MSE values help us see the model’s performance, while the coefficients reveal the impact of regularization.

Feel free to tweak the $( \alpha )$ values or add more features to the dataset to explore further!

Solving the Quadratic Assignment Problem

May 29, 2025

A Complete Guide with Python Implementation

The Quadratic Assignment Problem (QAP) is one of the most challenging combinatorial optimization problems in operations research.
Unlike simple assignment problems, QAP considers both the assignment costs and the interaction costs between facilities, making it particularly relevant for real-world facility location problems.

Problem Definition

The QAP can be mathematically formulated as:

$$\min \sum_{i=1}^{n} \sum_{j=1}^{n} \sum_{k=1}^{n} \sum_{l=1}^{n} f_{ik} \cdot d_{jl} \cdot x_{ij} \cdot x_{kl}$$

Subject to:
$$\sum_{j=1}^{n} x_{ij} = 1 \quad \forall i$$
$$\sum_{i=1}^{n} x_{ij} = 1 \quad \forall j$$
$$x_{ij} \in {0, 1} \quad \forall i,j$$

Where:

$f_{ik}$ is the flow between facilities $i$ and $k$
$d_{jl}$ is the distance between locations $j$ and $l$
$x_{ij} = 1$ if facility $i$ is assigned to location $j$, 0 otherwise

Practical Example: Hospital Department Layout

Let’s consider a hospital where we need to assign 4 departments to 4 available locations. The departments are:

Emergency Room (ER)
Surgery Department
Laboratory
Pharmacy

We need to minimize the total transportation cost considering both patient flow between departments and distances between locations.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from itertools import permutations
import pandas as pd
from scipy.optimize import linear_sum_assignment
import networkx as nx
import math

class QAPSolver:
    def __init__(self, flow_matrix, distance_matrix, facility_names, location_names):
        """
        Initialize QAP Solver
        
        Parameters:
        flow_matrix: Matrix representing flow between facilities
        distance_matrix: Matrix representing distances between locations  
        facility_names: Names of facilities
        location_names: Names of locations
        """
        self.flow_matrix = np.array(flow_matrix)
        self.distance_matrix = np.array(distance_matrix)
        self.facility_names = facility_names
        self.location_names = location_names
        self.n = len(facility_names)
        self.best_assignment = None
        self.best_cost = float('inf')
        self.all_solutions = []
        
    def calculate_cost(self, assignment):
        """
        Calculate total cost for a given assignment
        
        Parameters:
        assignment: List where assignment[i] is the location assigned to facility i
        
        Returns:
        Total cost of the assignment
        """
        total_cost = 0
        for i in range(self.n):
            for j in range(self.n):
                # Flow between facilities i and j
                flow = self.flow_matrix[i][j]
                # Distance between assigned locations
                distance = self.distance_matrix[assignment[i]][assignment[j]]
                total_cost += flow * distance
        return total_cost
    
    def solve_exhaustive(self):
        """
        Solve QAP using exhaustive search (suitable for small instances)
        """
        print("Solving QAP using exhaustive search...")
        
        # Generate all possible permutations
        for perm in permutations(range(self.n)):
            cost = self.calculate_cost(perm)
            self.all_solutions.append((list(perm), cost))
            
            if cost < self.best_cost:
                self.best_cost = cost
                self.best_assignment = list(perm)
        
        # Sort solutions by cost
        self.all_solutions.sort(key=lambda x: x[1])
        
        print(f"Optimal cost: {self.best_cost}")
        print(f"Optimal assignment: {self.best_assignment}")
        return self.best_assignment, self.best_cost
    
    def solve_greedy_heuristic(self):
        """
        Solve QAP using a greedy heuristic approach
        """
        print("Solving QAP using greedy heuristic...")
        
        # Calculate assignment priorities based on flow*distance products
        assignment_costs = np.zeros((self.n, self.n))
        
        for i in range(self.n):
            for j in range(self.n):
                # Calculate cost if facility i is assigned to location j
                cost = 0
                for k in range(self.n):
                    for l in range(self.n):
                        if k != i:  # Other facilities
                            # Estimate cost assuming worst case for other assignments  
                            flow = self.flow_matrix[i][k] + self.flow_matrix[k][i]
                            distance = self.distance_matrix[j][l]
                            cost += flow * distance
                assignment_costs[i][j] = cost
        
        # Use Hungarian algorithm for approximation
        row_ind, col_ind = linear_sum_assignment(assignment_costs)
        heuristic_assignment = col_ind.tolist()
        heuristic_cost = self.calculate_cost(heuristic_assignment)
        
        print(f"Heuristic cost: {heuristic_cost}")
        print(f"Heuristic assignment: {heuristic_assignment}")
        
        return heuristic_assignment, heuristic_cost
    
    def display_results(self):
        """Display the results in a formatted way"""
        print("\n" + "="*60)
        print("QUADRATIC ASSIGNMENT PROBLEM RESULTS")
        print("="*60)
        
        print(f"\nProblem Size: {self.n} facilities, {self.n} locations")
        print(f"Total possible assignments: {math.factorial(self.n)}")
        
        print(f"\nOptimal Assignment:")
        for i, loc in enumerate(self.best_assignment):
            print(f"  {self.facility_names[i]} → {self.location_names[loc]}")
        
        print(f"\nOptimal Cost: {self.best_cost}")
        
        # Show top 5 best and worst solutions
        print(f"\nTop 5 Best Solutions:")
        for i, (assignment, cost) in enumerate(self.all_solutions[:5]):
            assignment_str = " → ".join([f"{self.facility_names[j]}:{self.location_names[assignment[j]]}" 
                                       for j in range(self.n)])
            print(f"  {i+1}. Cost: {cost:6.1f} | {assignment_str}")
            
        print(f"\nTop 5 Worst Solutions:")
        for i, (assignment, cost) in enumerate(self.all_solutions[-5:]):
            assignment_str = " → ".join([f"{self.facility_names[j]}:{self.location_names[assignment[j]]}" 
                                       for j in range(self.n)])
            print(f"  {i+1}. Cost: {cost:6.1f} | {assignment_str}")

    def visualize_results(self):
        """Create comprehensive visualizations of the QAP solution"""
        
        # Set up the plotting style
        plt.style.use('default')
        sns.set_palette("husl")
        
        # Create figure with subplots
        fig = plt.figure(figsize=(20, 15))
        
        # 1. Flow Matrix Heatmap
        ax1 = plt.subplot(2, 4, 1)
        sns.heatmap(self.flow_matrix, annot=True, fmt='g', cmap='YlOrRd',
                   xticklabels=self.facility_names, yticklabels=self.facility_names,
                   cbar_kws={'label': 'Flow Units'})
        plt.title('Flow Matrix Between Facilities', fontsize=12, fontweight='bold')
        plt.xlabel('To Facility')
        plt.ylabel('From Facility')
        
        # 2. Distance Matrix Heatmap  
        ax2 = plt.subplot(2, 4, 2)
        sns.heatmap(self.distance_matrix, annot=True, fmt='g', cmap='Blues',
                   xticklabels=self.location_names, yticklabels=self.location_names,
                   cbar_kws={'label': 'Distance Units'})
        plt.title('Distance Matrix Between Locations', fontsize=12, fontweight='bold')
        plt.xlabel('To Location')
        plt.ylabel('From Location')
        
        # 3. Cost Distribution
        ax3 = plt.subplot(2, 4, 3)
        costs = [sol[1] for sol in self.all_solutions]
        plt.hist(costs, bins=min(20, len(costs)//2), alpha=0.7, color='skyblue', edgecolor='black')
        plt.axvline(self.best_cost, color='red', linestyle='--', linewidth=2, label=f'Optimal: {self.best_cost}')
        plt.xlabel('Total Cost')
        plt.ylabel('Number of Solutions')
        plt.title('Distribution of Solution Costs', fontsize=12, fontweight='bold')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # 4. Assignment Network Graph
        ax4 = plt.subplot(2, 4, 4)
        G = nx.Graph()
        
        # Add nodes for facilities and locations
        facility_nodes = [f"F_{i}_{self.facility_names[i][:3]}" for i in range(self.n)]
        location_nodes = [f"L_{i}_{self.location_names[i][:3]}" for i in range(self.n)]
        
        G.add_nodes_from(facility_nodes)
        G.add_nodes_from(location_nodes)
        
        # Add edges for optimal assignment
        for i, loc in enumerate(self.best_assignment):
            G.add_edge(facility_nodes[i], location_nodes[loc])
        
        # Position nodes
        pos = {}
        for i, node in enumerate(facility_nodes):
            pos[node] = (0, i)
        for i, node in enumerate(location_nodes):
            pos[node] = (2, i)
            
        nx.draw(G, pos, ax=ax4, with_labels=True, node_color=['lightcoral']*self.n + ['lightblue']*self.n,
                node_size=1500, font_size=8, font_weight='bold')
        plt.title('Optimal Assignment Network', fontsize=12, fontweight='bold')
        
        # 5. Cost Contribution Matrix
        ax5 = plt.subplot(2, 4, 5)
        cost_matrix = np.zeros((self.n, self.n))
        for i in range(self.n):
            for j in range(self.n):
                flow = self.flow_matrix[i][j]
                distance = self.distance_matrix[self.best_assignment[i]][self.best_assignment[j]]
                cost_matrix[i][j] = flow * distance
                
        sns.heatmap(cost_matrix, annot=True, fmt='.1f', cmap='Reds',
                   xticklabels=self.facility_names, yticklabels=self.facility_names,
                   cbar_kws={'label': 'Cost Contribution'})
        plt.title('Cost Contribution Matrix\n(Optimal Assignment)', fontsize=12, fontweight='bold')
        plt.xlabel('To Facility')
        plt.ylabel('From Facility')
        
        # 6. Solution Quality Comparison
        ax6 = plt.subplot(2, 4, 6)
        x = range(len(self.all_solutions))
        y = [sol[1] for sol in self.all_solutions]
        plt.plot(x, y, 'b-', alpha=0.7, linewidth=1, label='All Solutions')
        plt.scatter(0, self.best_cost, color='red', s=100, zorder=5, label=f'Optimal: {self.best_cost}')
        plt.xlabel('Solution Rank')
        plt.ylabel('Total Cost')
        plt.title('Solution Quality Ranking', fontsize=12, fontweight='bold')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # 7. Facility-Location Assignment Bar Chart
        ax7 = plt.subplot(2, 4, 7)
        y_pos = np.arange(self.n)
        colors = sns.color_palette("husl", self.n)
        
        bars = plt.barh(y_pos, [self.best_assignment[i] for i in range(self.n)], color=colors)
        plt.yticks(y_pos, self.facility_names)
        plt.xlabel('Assigned Location Index')
        plt.title('Optimal Facility Assignments', fontsize=12, fontweight='bold')
        plt.grid(True, alpha=0.3, axis='x')
        
        # Add location names as text
        for i, bar in enumerate(bars):
            width = bar.get_width()
            location_name = self.location_names[int(width)]
            plt.text(width + 0.05, bar.get_y() + bar.get_height()/2, 
                    location_name, ha='left', va='center', fontweight='bold')
        
        # 8. Performance Metrics
        ax8 = plt.subplot(2, 4, 8)
        ax8.axis('off')
        
        # Calculate some performance metrics
        total_flow = np.sum(self.flow_matrix)
        total_distance = np.sum(self.distance_matrix)
        avg_cost = np.mean(costs)
        std_cost = np.std(costs)
        worst_cost = max(costs)
        gap = ((worst_cost - self.best_cost) / self.best_cost) * 100
        
        metrics_text = f"""
        PERFORMANCE METRICS
        ─────────────────────
        Total Flow: {total_flow:.1f}
        Total Distance: {total_distance:.1f}
        
        Optimal Cost: {self.best_cost:.1f}
        Average Cost: {avg_cost:.1f}
        Worst Cost: {worst_cost:.1f}
        Standard Deviation: {std_cost:.1f}
        
        Optimality Gap: {gap:.1f}%
        Solutions Evaluated: {len(self.all_solutions)}
        """
        
        plt.text(0.1, 0.9, metrics_text, transform=ax8.transAxes, fontsize=11,
                verticalalignment='top', fontfamily='monospace',
                bbox=dict(boxstyle="round,pad=0.3", facecolor="lightgray", alpha=0.8))
        
        plt.tight_layout()
        plt.suptitle('Quadratic Assignment Problem: Complete Analysis', 
                    fontsize=16, fontweight='bold', y=0.98)
        plt.show()
        
        # Create a detailed assignment comparison table
        self.create_assignment_table()
    
    def create_assignment_table(self):
        """Create a detailed table comparing different assignments"""
        print("\n" + "="*80)
        print("DETAILED ASSIGNMENT ANALYSIS")
        print("="*80)
        
        # Create DataFrame for analysis
        analysis_data = []
        
        for rank, (assignment, cost) in enumerate(self.all_solutions[:10]):  # Top 10 solutions
            row = {
                'Rank': rank + 1,
                'Total_Cost': cost,
                'Cost_Gap_%': ((cost - self.best_cost) / self.best_cost) * 100
            }
            
            # Add individual assignments
            for i, loc in enumerate(assignment):
                row[f'{self.facility_names[i]}'] = self.location_names[loc]
            
            analysis_data.append(row)
            
        df = pd.DataFrame(analysis_data)
        print(df.to_string(index=False))
        
        # Flow-Distance Analysis
        print(f"\n\nFLOW-DISTANCE INTERACTION ANALYSIS (Optimal Solution)")
        print("-" * 60)
        
        for i in range(self.n):
            for j in range(self.n):
                if i != j:  # Skip diagonal
                    flow = self.flow_matrix[i][j]
                    if flow > 0:  # Only show non-zero flows
                        loc_i = self.best_assignment[i]
                        loc_j = self.best_assignment[j]
                        distance = self.distance_matrix[loc_i][loc_j]
                        cost_contribution = flow * distance
                        
                        print(f"{self.facility_names[i]:>12} → {self.facility_names[j]:<12} | "
                              f"Flow: {flow:>4.1f} × Distance: {distance:>4.1f} = Cost: {cost_contribution:>6.1f}")

# Define the hospital department layout problem
def create_hospital_qap():
    """
    Create a hospital department layout QAP instance
    """
    
    # Flow matrix (patients/staff movement per day between departments)
    # Rows/Cols: ER, Surgery, Laboratory, Pharmacy
    flow_matrix = [
        [0,   25,  30,  40],  # ER to others
        [25,  0,   35,  20],  # Surgery to others  
        [30,  35,  0,   15],  # Laboratory to others
        [40,  20,  15,  0]    # Pharmacy to others
    ]
    
    # Distance matrix between locations (in meters)
    # Rows/Cols: Location A, B, C, D
    distance_matrix = [
        [0,   50,  80,  120], # Location A to others
        [50,  0,   60,  90],  # Location B to others
        [80,  60,  0,   40],  # Location C to others  
        [120, 90,  40,  0]    # Location D to others
    ]
    
    facility_names = ['Emergency Room', 'Surgery Dept', 'Laboratory', 'Pharmacy']
    location_names = ['Location A', 'Location B', 'Location C', 'Location D'] 
    
    return QAPSolver(flow_matrix, distance_matrix, facility_names, location_names)

# Main execution
if __name__ == "__main__":
    print("QUADRATIC ASSIGNMENT PROBLEM SOLVER")
    print("Hospital Department Layout Optimization")
    print("="*50)
    
    # Create and solve the hospital QAP instance
    qap = create_hospital_qap()
    
    # Solve using exhaustive search
    optimal_assignment, optimal_cost = qap.solve_exhaustive()
    
    # Solve using heuristic for comparison  
    heuristic_assignment, heuristic_cost = qap.solve_greedy_heuristic()
    
    # Display results
    qap.display_results()
    
    # Compare heuristic vs optimal
    if optimal_cost != heuristic_cost:
        gap = ((heuristic_cost - optimal_cost) / optimal_cost) * 100
        print(f"\nHeuristic vs Optimal Comparison:")
        print(f"Optimal Cost: {optimal_cost}")
        print(f"Heuristic Cost: {heuristic_cost}")
        print(f"Gap: {gap:.2f}%")
    else:
        print(f"\nHeuristic found the optimal solution!")
    
    # Create visualizations
    qap.visualize_results()

Detailed Code Explanation

Class Structure and Initialization

The QAPSolver class is designed to handle quadratic assignment problems efficiently.
The constructor takes four key parameters:

Flow Matrix: Represents the interaction intensity between facilities (e.g., patient movement between hospital departments)
Distance Matrix: Represents the physical distances between available locations
Facility Names: Human-readable names for better interpretation
Location Names: Human-readable location identifiers

Core Algorithm Implementation

1. Cost Calculation Function

The calculate_cost() method implements the QAP objective function:

def calculate_cost(self, assignment):
    total_cost = 0
    for i in range(self.n):
        for j in range(self.n):
            flow = self.flow_matrix[i][j]
            distance = self.distance_matrix[assignment[i]][assignment[j]]
            total_cost += flow * distance
    return total_cost

This function computes the total cost by summing all flow-distance products for every pair of facilities in their assigned locations.

2. Exhaustive Search Solution

The solve_exhaustive() method generates all possible permutations of facility assignments.
For our 4×4 problem, this means evaluating 4! = 24 different assignments.
While computationally expensive for larger instances, this guarantees finding the global optimum for small problems.

3. Greedy Heuristic Approach

The solve_greedy_heuristic() method provides a fast approximation using the Hungarian algorithm.
It creates an assignment cost matrix and solves it as a linear assignment problem, which often provides good solutions quickly.

Visualization System

The visualization system creates eight different plots to provide comprehensive insights:

Flow Matrix Heatmap: Shows interaction intensities between facilities
Distance Matrix Heatmap: Displays physical distances between locations
Cost Distribution: Histogram showing the distribution of all solution costs
Assignment Network: Graph representation of the optimal assignment
Cost Contribution Matrix: Shows how much each facility pair contributes to total cost
Solution Quality Ranking: Line plot showing all solutions ranked by cost
Facility Assignment Bar Chart: Visual representation of which facility goes where
Performance Metrics: Summary statistics and key performance indicators

Expected Results and Interpretation

When you run this code, you’ll see:

QUADRATIC ASSIGNMENT PROBLEM SOLVER
Hospital Department Layout Optimization
==================================================
Solving QAP using exhaustive search...
Optimal cost: 21700
Optimal assignment: [2, 1, 0, 3]
Solving QAP using greedy heuristic...
Heuristic cost: 21700
Heuristic assignment: [2, 1, 0, 3]

============================================================
QUADRATIC ASSIGNMENT PROBLEM RESULTS
============================================================

Problem Size: 4 facilities, 4 locations
Total possible assignments: 24

Optimal Assignment:
  Emergency Room → Location C
  Surgery Dept → Location B
  Laboratory → Location A
  Pharmacy → Location D

Optimal Cost: 21700

Top 5 Best Solutions:
  1. Cost: 21700.0 | Emergency Room:Location C → Surgery Dept:Location B → Laboratory:Location A → Pharmacy:Location D
  2. Cost: 21800.0 | Emergency Room:Location C → Surgery Dept:Location A → Laboratory:Location B → Pharmacy:Location D
  3. Cost: 22000.0 | Emergency Room:Location B → Surgery Dept:Location C → Laboratory:Location D → Pharmacy:Location A
  4. Cost: 22100.0 | Emergency Room:Location B → Surgery Dept:Location D → Laboratory:Location C → Pharmacy:Location A
  5. Cost: 23000.0 | Emergency Room:Location A → Surgery Dept:Location D → Laboratory:Location C → Pharmacy:Location B

Top 5 Worst Solutions:
  1. Cost: 25500.0 | Emergency Room:Location B → Surgery Dept:Location A → Laboratory:Location D → Pharmacy:Location C
  2. Cost: 25900.0 | Emergency Room:Location A → Surgery Dept:Location B → Laboratory:Location C → Pharmacy:Location D
  3. Cost: 25900.0 | Emergency Room:Location D → Surgery Dept:Location C → Laboratory:Location A → Pharmacy:Location B
  4. Cost: 25900.0 | Emergency Room:Location D → Surgery Dept:Location C → Laboratory:Location B → Pharmacy:Location A
  5. Cost: 26000.0 | Emergency Room:Location A → Surgery Dept:Location B → Laboratory:Location D → Pharmacy:Location C

Heuristic found the optimal solution!


================================================================================
DETAILED ASSIGNMENT ANALYSIS
================================================================================
 Rank  Total_Cost  Cost_Gap_% Emergency Room Surgery Dept Laboratory   Pharmacy
    1       21700    0.000000     Location C   Location B Location A Location D
    2       21800    0.460829     Location C   Location A Location B Location D
    3       22000    1.382488     Location B   Location C Location D Location A
    4       22100    1.843318     Location B   Location D Location C Location A
    5       23000    5.990783     Location A   Location D Location C Location B
    6       23100    6.451613     Location A   Location C Location D Location B
    7       23100    6.451613     Location D   Location A Location B Location C
    8       23200    6.912442     Location D   Location B Location A Location C
    9       23700    9.216590     Location C   Location B Location D Location A
   10       24000   10.599078     Location B   Location C Location A Location D


FLOW-DISTANCE INTERACTION ANALYSIS (Optimal Solution)
------------------------------------------------------------
Emergency Room → Surgery Dept | Flow: 25.0 × Distance: 60.0 = Cost: 1500.0
Emergency Room → Laboratory   | Flow: 30.0 × Distance: 80.0 = Cost: 2400.0
Emergency Room → Pharmacy     | Flow: 40.0 × Distance: 40.0 = Cost: 1600.0
Surgery Dept → Emergency Room | Flow: 25.0 × Distance: 60.0 = Cost: 1500.0
Surgery Dept → Laboratory   | Flow: 35.0 × Distance: 50.0 = Cost: 1750.0
Surgery Dept → Pharmacy     | Flow: 20.0 × Distance: 90.0 = Cost: 1800.0
  Laboratory → Emergency Room | Flow: 30.0 × Distance: 80.0 = Cost: 2400.0
  Laboratory → Surgery Dept | Flow: 35.0 × Distance: 50.0 = Cost: 1750.0
  Laboratory → Pharmacy     | Flow: 15.0 × Distance: 120.0 = Cost: 1800.0
    Pharmacy → Emergency Room | Flow: 40.0 × Distance: 40.0 = Cost: 1600.0
    Pharmacy → Surgery Dept | Flow: 20.0 × Distance: 90.0 = Cost: 1800.0
    Pharmacy → Laboratory   | Flow: 15.0 × Distance: 120.0 = Cost: 1800.0

Optimal Solution Analysis

The algorithm will find the assignment that minimizes total transportation cost.
For our hospital example, this typically involves:

Placing high-flow facility pairs (like ER-Pharmacy) close together
Minimizing the sum of flow×distance products across all facility pairs

Performance Metrics

The code provides several key metrics:

Optimality Gap: Difference between best and worst solutions
Solution Distribution: How costs are distributed across all possible assignments
Cost Contributions: Which facility pairs contribute most to total cost

Visual Insights

The comprehensive visualization reveals:

Hot Spots: High-flow corridors in the optimal layout
Distance Efficiency: How well the solution utilizes short distances for high flows
Trade-offs: Where the algorithm chose to place lower-priority flows at longer distances

Mathematical Foundation

The QAP’s complexity comes from its quadratic nature - the cost depends on products of two decision variables ($x_{ij} \cdot x_{kl}$). This makes it:

NP-hard: No polynomial-time algorithm is known
Non-linear: Cannot be solved with standard linear programming
Combinatorially explosive: n! possible solutions for n facilities

Practical Applications

This QAP formulation applies to many real-world scenarios:

Manufacturing: Assembly line layout optimization
Computing: Processor task scheduling
Urban Planning: Public service facility placement
Supply Chain: Warehouse and distribution center location

The hospital example demonstrates how QAP captures both direct assignment costs and interaction effects, making it particularly valuable for facility layout problems where workflows and distances both matter significantly.

Demand Forecasting and Inventory Optimization

May 28, 2025

A Machine Learning + Operations Research Approach

Today we’re diving into one of the most fascinating applications where machine learning meets operations research: demand forecasting combined with inventory optimization.
We’ll tackle a real-world scenario where a retail company needs to predict future demand and optimize their inventory levels simultaneously.

The Problem Setup

Imagine you’re managing inventory for an electronics retailer that sells smartphones. You need to:

Predict future demand using historical sales data
Optimize inventory levels to minimize total costs (holding + stockout costs)
Account for uncertainty in demand predictions

The mathematical formulation combines:

ML Component: Time series forecasting with uncertainty quantification
OR Component: Stochastic inventory optimization

The objective function we want to minimize is:

$$\min_{Q_t} \mathbb{E}\left[\sum_{t=1}^{T} \left(h \cdot I_t^+ + p \cdot I_t^- + c \cdot Q_t\right)\right]$$

Where:

$Q_t$ = order quantity at time $t$
$I_t^+$ = positive inventory (holding cost)
$I_t^-$ = negative inventory (stockout cost)
$h$ = holding cost per unit
$p$ = penalty cost per unit shortage
$c$ = ordering cost per unit

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.preprocessing import StandardScaler
from scipy import optimize
import warnings
warnings.filterwarnings('ignore')

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

class DemandForecaster:
    """Machine Learning based demand forecaster with uncertainty quantification"""
    
    def __init__(self):
        self.model = GradientBoostingRegressor(
            n_estimators=100,
            learning_rate=0.1,
            max_depth=4,
            random_state=42
        )
        self.scaler = StandardScaler()
        self.is_fitted = False
        
    def create_features(self, data, lookback=7):
        """Create time series features for ML model"""
        features = []
        targets = []
        
        for i in range(lookback, len(data)):
            # Lagged demand features
            lag_features = data[i-lookback:i].tolist()
            
            # Additional features
            day_of_week = i % 7
            week_of_year = (i // 7) % 52
            trend = i / len(data)  # Linear trend
            
            # Seasonal features
            seasonal = np.sin(2 * np.pi * i / 7)  # Weekly seasonality
            
            feature_vector = lag_features + [day_of_week, week_of_year, trend, seasonal]
            features.append(feature_vector)
            targets.append(data[i])
            
        return np.array(features), np.array(targets)
    
    def fit(self, demand_data):
        """Train the demand forecasting model"""
        X, y = self.create_features(demand_data)
        X_scaled = self.scaler.fit_transform(X)
        self.model.fit(X_scaled, y)
        self.is_fitted = True
        return self
    
    def predict_with_uncertainty(self, demand_data, horizon=30, n_simulations=100):
        """Predict demand with uncertainty bounds using quantile regression approach"""
        if not self.is_fitted:
            raise ValueError("Model must be fitted before prediction")
            
        X, _ = self.create_features(demand_data)
        X_scaled = self.scaler.transform(X[-1:])  # Use last observation
        
        predictions = []
        uncertainties = []
        
        # Generate probabilistic forecasts
        for h in range(horizon):
            # Base prediction
            base_pred = self.model.predict(X_scaled)[0]
            
            # Add uncertainty based on model's training error patterns
            # In practice, you'd use more sophisticated uncertainty quantification
            residual_std = np.std(demand_data) * (1 + 0.1 * h)  # Increasing uncertainty over time
            pred_samples = np.random.normal(base_pred, residual_std, n_simulations)
            pred_samples = np.maximum(pred_samples, 0)  # Demand can't be negative
            
            predictions.append({
                'mean': base_pred,
                'median': np.median(pred_samples),
                'lower_80': np.percentile(pred_samples, 10),
                'upper_80': np.percentile(pred_samples, 90),
                'lower_95': np.percentile(pred_samples, 2.5),
                'upper_95': np.percentile(pred_samples, 97.5),
                'samples': pred_samples
            })
            
            # Update features for next prediction (simplified)
            # In practice, you'd properly roll the window
            X_scaled = np.roll(X_scaled, -1)
            X_scaled[0, -4:] = [h % 7, (h // 7) % 52, h / horizon, np.sin(2 * np.pi * h / 7)]
        
        return predictions

class InventoryOptimizer:
    """Operations Research based inventory optimizer"""
    
    def __init__(self, holding_cost=1.0, stockout_cost=10.0, order_cost=0.5):
        self.holding_cost = holding_cost
        self.stockout_cost = stockout_cost
        self.order_cost = order_cost
    
    def newsvendor_optimal_quantity(self, demand_samples):
        """Calculate optimal order quantity using newsvendor model"""
        # Critical ratio for newsvendor problem
        critical_ratio = self.stockout_cost / (self.stockout_cost + self.holding_cost)
        optimal_quantity = np.percentile(demand_samples, critical_ratio * 100)
        return max(0, optimal_quantity)
    
    def calculate_expected_cost(self, order_quantity, demand_samples):
        """Calculate expected total cost for given order quantity"""
        inventory_levels = order_quantity - demand_samples
        
        # Holding costs (positive inventory)
        holding_costs = np.maximum(inventory_levels, 0) * self.holding_cost
        
        # Stockout costs (negative inventory)
        stockout_costs = np.maximum(-inventory_levels, 0) * self.stockout_cost
        
        # Order costs
        order_costs = order_quantity * self.order_cost
        
        total_cost = np.mean(holding_costs + stockout_costs) + order_costs
        return total_cost
    
    def optimize_inventory(self, demand_predictions, initial_inventory=0):
        """Optimize inventory decisions over planning horizon"""
        horizon = len(demand_predictions)
        optimal_orders = []
        expected_costs = []
        inventory_levels = [initial_inventory]
        
        current_inventory = initial_inventory
        
        for t in range(horizon):
            demand_samples = demand_predictions[t]['samples']
            
            # Calculate net demand (demand - current inventory)
            net_demand_samples = np.maximum(demand_samples - current_inventory, 0)
            
            # Find optimal order quantity for this period
            if len(net_demand_samples[net_demand_samples > 0]) > 0:
                optimal_qty = self.newsvendor_optimal_quantity(net_demand_samples)
            else:
                optimal_qty = 0
                
            optimal_orders.append(optimal_qty)
            
            # Calculate expected cost
            total_inventory = current_inventory + optimal_qty
            cost = self.calculate_expected_cost(optimal_qty, demand_samples)
            expected_costs.append(cost)
            
            # Update inventory for next period (simplified simulation)
            expected_demand = demand_predictions[t]['mean']
            current_inventory = max(0, total_inventory - expected_demand)
            inventory_levels.append(current_inventory)
        
        return {
            'optimal_orders': optimal_orders,
            'expected_costs': expected_costs,
            'inventory_levels': inventory_levels[:-1]  # Remove last element
        }

# Generate synthetic demand data
def generate_demand_data(n_periods=200):
    """Generate realistic demand data with trend, seasonality, and randomness"""
    np.random.seed(42)
    
    # Base demand level
    base_demand = 100
    
    # Trend component
    trend = np.linspace(0, 20, n_periods)
    
    # Seasonal components
    weekly_season = 30 * np.sin(2 * np.pi * np.arange(n_periods) / 7)
    monthly_season = 15 * np.sin(2 * np.pi * np.arange(n_periods) / 30)
    
    # Random noise
    noise = np.random.normal(0, 10, n_periods)
    
    # Combine components
    demand = base_demand + trend + weekly_season + monthly_season + noise
    demand = np.maximum(demand, 0)  # Ensure non-negative demand
    
    return demand

# Main execution
print("🚀 Starting Integrated Demand Forecasting and Inventory Optimization")
print("=" * 70)

# Step 1: Generate and prepare data
demand_data = generate_demand_data(150)
train_size = int(0.8 * len(demand_data))
train_data = demand_data[:train_size]
test_data = demand_data[train_size:]

print(f"📊 Data Summary:")
print(f"   Total periods: {len(demand_data)}")
print(f"   Training periods: {len(train_data)}")
print(f"   Test periods: {len(test_data)}")
print(f"   Average demand: {np.mean(demand_data):.2f}")
print(f"   Demand std: {np.std(demand_data):.2f}")

# Step 2: Train demand forecasting model
print("\n🤖 Training Machine Learning Demand Forecaster...")
forecaster = DemandForecaster()
forecaster.fit(train_data)

# Step 3: Generate demand predictions
print("🔮 Generating Demand Predictions with Uncertainty...")
forecast_horizon = 30
demand_predictions = forecaster.predict_with_uncertainty(
    train_data, 
    horizon=forecast_horizon,
    n_simulations=1000
)

# Step 4: Optimize inventory
print("⚙️ Optimizing Inventory Decisions...")
optimizer = InventoryOptimizer(
    holding_cost=2.0,     # $2 per unit per period
    stockout_cost=15.0,   # $15 penalty per unit shortage
    order_cost=1.0        # $1 per unit ordered
)

optimization_results = optimizer.optimize_inventory(demand_predictions)

# Step 5: Performance Analysis
print("\n📈 Performance Analysis:")
print("-" * 30)

# Forecasting accuracy (using available test data)
test_predictions = []
for i in range(min(len(test_data), forecast_horizon)):
    test_predictions.append(demand_predictions[i]['mean'])

if len(test_predictions) > 0:
    mae = mean_absolute_error(test_data[:len(test_predictions)], test_predictions)
    rmse = np.sqrt(mean_squared_error(test_data[:len(test_predictions)], test_predictions))
    mape = np.mean(np.abs((test_data[:len(test_predictions)] - test_predictions) / test_data[:len(test_predictions)])) * 100
    
    print(f"Forecasting Performance:")
    print(f"  MAE: {mae:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAPE: {mape:.2f}%")

# Inventory optimization results
total_expected_cost = sum(optimization_results['expected_costs'])
avg_order_quantity = np.mean(optimization_results['optimal_orders'])
avg_inventory_level = np.mean(optimization_results['inventory_levels'])

print(f"\nInventory Optimization Results:")
print(f"  Total expected cost: ${total_expected_cost:.2f}")
print(f"  Average order quantity: {avg_order_quantity:.2f} units")
print(f"  Average inventory level: {avg_inventory_level:.2f} units")

print("\n✅ Analysis Complete! Generating visualizations...")

# Create comprehensive visualizations
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Integrated Demand Forecasting and Inventory Optimization Results', fontsize=16, fontweight='bold')

# Plot 1: Historical demand and forecasts
ax1 = axes[0, 0]
time_historical = range(len(demand_data))
time_forecast = range(len(train_data), len(train_data) + forecast_horizon)

ax1.plot(time_historical, demand_data, 'b-', alpha=0.7, label='Historical Demand', linewidth=2)
ax1.axvline(x=len(train_data), color='red', linestyle='--', alpha=0.7, label='Train/Test Split')

# Plot forecast with uncertainty bands
forecast_means = [pred['mean'] for pred in demand_predictions]
forecast_lower_80 = [pred['lower_80'] for pred in demand_predictions]
forecast_upper_80 = [pred['upper_80'] for pred in demand_predictions]
forecast_lower_95 = [pred['lower_95'] for pred in demand_predictions]
forecast_upper_95 = [pred['upper_95'] for pred in demand_predictions]

ax1.plot(time_forecast, forecast_means, 'g-', linewidth=2, label='Forecast Mean')
ax1.fill_between(time_forecast, forecast_lower_80, forecast_upper_80, alpha=0.3, color='green', label='80% Confidence')
ax1.fill_between(time_forecast, forecast_lower_95, forecast_upper_95, alpha=0.1, color='green', label='95% Confidence')

ax1.set_title('Demand Forecasting with Uncertainty Quantification')
ax1.set_xlabel('Time Period')
ax1.set_ylabel('Demand')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Optimal order quantities and inventory levels
ax2 = axes[0, 1]
forecast_periods = range(forecast_horizon)

ax2_twin = ax2.twinx()
bars = ax2.bar(forecast_periods, optimization_results['optimal_orders'], alpha=0.6, color='orange', label='Optimal Orders')
line = ax2_twin.plot(forecast_periods, optimization_results['inventory_levels'], 'r-', linewidth=2, marker='o', label='Inventory Level')

ax2.set_title('Optimal Ordering Policy and Inventory Levels')
ax2.set_xlabel('Time Period')
ax2.set_ylabel('Order Quantity', color='orange')
ax2_twin.set_ylabel('Inventory Level', color='red')
ax2.grid(True, alpha=0.3)

# Create combined legend
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, loc='upper right')

# Plot 3: Cost breakdown analysis
ax3 = axes[1, 0]
cost_components = []
periods = range(forecast_horizon)

for i, pred in enumerate(demand_predictions):
    order_qty = optimization_results['optimal_orders'][i]
    demand_samples = pred['samples']
    
    # Calculate cost components
    inventory_after_order = optimization_results['inventory_levels'][i] + order_qty
    inventory_levels = inventory_after_order - demand_samples
    
    holding_cost = np.mean(np.maximum(inventory_levels, 0)) * optimizer.holding_cost
    stockout_cost = np.mean(np.maximum(-inventory_levels, 0)) * optimizer.stockout_cost
    order_cost = order_qty * optimizer.order_cost
    
    cost_components.append([holding_cost, stockout_cost, order_cost])

cost_components = np.array(cost_components)

ax3.stackplot(periods, 
              cost_components[:, 0], 
              cost_components[:, 1], 
              cost_components[:, 2],
              labels=['Holding Cost', 'Stockout Cost', 'Order Cost'],
              alpha=0.7)

ax3.set_title('Cost Component Analysis Over Time')
ax3.set_xlabel('Time Period')
ax3.set_ylabel('Expected Cost ($)')
ax3.legend(loc='upper left')
ax3.grid(True, alpha=0.3)

# Plot 4: Demand distribution and optimal policy
ax4 = axes[1, 1]
# Show distribution for first few periods
sample_period = 5
demand_samples = demand_predictions[sample_period]['samples']
optimal_order = optimization_results['optimal_orders'][sample_period]

ax4.hist(demand_samples, bins=30, alpha=0.7, density=True, color='skyblue', label='Demand Distribution')
ax4.axvline(x=np.mean(demand_samples), color='blue', linestyle='-', linewidth=2, label=f'Mean Demand: {np.mean(demand_samples):.1f}')
ax4.axvline(x=optimal_order, color='red', linestyle='--', linewidth=2, label=f'Optimal Order: {optimal_order:.1f}')

ax4.set_title(f'Demand Distribution and Optimal Policy (Period {sample_period+1})')
ax4.set_xlabel('Demand')
ax4.set_ylabel('Probability Density')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Service level and fill rate
print("\n📊 Additional Performance Metrics:")
print("-" * 40)

service_levels = []
fill_rates = []

for i, pred in enumerate(demand_predictions):
    order_qty = optimization_results['optimal_orders'][i]
    current_inv = optimization_results['inventory_levels'][i]
    total_available = current_inv + order_qty
    demand_samples = pred['samples']
    
    # Service level (probability of no stockout)
    service_level = np.mean(demand_samples <= total_available)
    service_levels.append(service_level)
    
    # Fill rate (expected fraction of demand satisfied)
    satisfied_demand = np.minimum(demand_samples, total_available)
    fill_rate = np.mean(satisfied_demand / demand_samples)
    fill_rates.append(fill_rate)

avg_service_level = np.mean(service_levels)
avg_fill_rate = np.mean(fill_rates)

print(f"Average Service Level: {avg_service_level:.1%}")
print(f"Average Fill Rate: {avg_fill_rate:.1%}")
print(f"Total Inventory Investment: ${np.sum(optimization_results['optimal_orders']):.2f}")

# Summary insights
print("\n🎯 Key Insights:")
print("-" * 20)
print("• The ML model captures both trend and seasonal patterns in demand")
print("• Uncertainty quantification enables robust inventory decisions")
print("• The newsvendor model balances holding and stockout costs optimally")
print(f"• Service level of {avg_service_level:.1%} achieved with optimized policy")
print("• Integration of ML and OR provides superior performance than using either alone")

print("\n🏁 Analysis completed successfully!")

Code Deep Dive

Let me break down the key components of this integrated solution:

1. The DemandForecaster Class

This class implements our machine learning component using gradient boosting:

1	def create_features(self, data, lookback=7):

The feature engineering creates:

Lagged features: Past 7 days of demand (autoregressive terms)
Calendar features: Day of week, week of year
Trend features: Linear time trend
Seasonal features: Weekly seasonality using sine waves

The mathematical representation of our features is:

2. Uncertainty Quantification

The predict_with_uncertainty method is crucial:

1 2	residual_std = np.std(demand_data) * (1 + 0.1 * h) pred_samples = np.random.normal(base_pred, residual_std, n_simulations)

This generates probabilistic forecasts by:

Creating multiple demand scenarios
Increasing uncertainty over longer horizons
Providing confidence intervals for decision-making

3. The InventoryOptimizer Class

This implements the operations research component:

The newsvendor model finds the optimal order quantity $Q^*$ that satisfies:
$$P(D \leq Q^*) = \frac{p}{p + h}$$

Where the critical ratio balances stockout penalty $p$ and holding cost $h$.

1
2
3

def newsvendor_optimal_quantity(self, demand_samples):
    critical_ratio = self.stockout_cost / (self.stockout_cost + self.holding_cost)
    optimal_quantity = np.percentile(demand_samples, critical_ratio * 100)

4. Cost Function Implementation

The expected total cost calculation:
$$E[C(Q)] = h \cdot E[\max(Q-D, 0)] + p \cdot E[\max(D-Q, 0)] + c \cdot Q$$

def calculate_expected_cost(self, order_quantity, demand_samples):
    inventory_levels = order_quantity - demand_samples
    holding_costs = np.maximum(inventory_levels, 0) * self.holding_cost
    stockout_costs = np.maximum(-inventory_levels, 0) * self.stockout_cost
    order_costs = order_quantity * self.order_cost

Results

🚀 Starting Integrated Demand Forecasting and Inventory Optimization
======================================================================
📊 Data Summary:
   Total periods: 150
   Training periods: 120
   Test periods: 30
   Average demand: 109.53
   Demand std: 25.80

🤖 Training Machine Learning Demand Forecaster...
🔮 Generating Demand Predictions with Uncertainty...
⚙️ Optimizing Inventory Decisions...

📈 Performance Analysis:
------------------------------
Forecasting Performance:
  MAE: 27.89
  RMSE: 30.94
  MAPE: 24.67%

Inventory Optimization Results:
  Total expected cost: $15434.82
  Average order quantity: 119.30 units
  Average inventory level: 69.18 units

✅ Analysis Complete! Generating visualizations...

Results Analysis

The visualization shows four key insights:

Plot 1: Forecasting Performance

The ML model captures both trend and weekly seasonality
Uncertainty bands widen over time (realistic prediction intervals)
The 80% and 95% confidence intervals provide decision-makers with risk quantification

Plot 2: Optimal Policy

Order quantities vary based on predicted demand and uncertainty
Inventory levels are maintained at service-level appropriate levels
The policy adapts to demand patterns dynamically

Plot 3: Cost Breakdown

Holding costs dominate when demand is low
Stockout costs spike during high-demand periods
Order costs remain relatively stable
This shows the trade-off optimization in action

Plot 4: Decision Illustration

Shows how the optimal order quantity positions itself relative to the demand distribution
The red dashed line (optimal order) is positioned based on the critical ratio
Higher stockout penalties push the optimal order to the right

Why This Integration Works

The key advantage of combining ML and OR is:

ML provides realistic demand scenarios with proper uncertainty quantification
OR provides optimal decisions given those scenarios
Together they handle both prediction and prescription

The mathematical beauty lies in the critical ratio formula:
$$\text{Service Level} = \frac{\text{Stockout Cost}}{\text{Stockout Cost + Holding Cost}}$$

This automatically balances risks based on business costs, not arbitrary service level targets.

Performance Metrics

The solution achieves:

Forecasting accuracy: MAE, RMSE, and MAPE metrics
Service level: Probability of avoiding stockouts
Fill rate: Percentage of demand satisfied
Cost optimization: Minimized total expected costs

This integrated approach outperforms traditional methods that treat forecasting and optimization separately, as it properly accounts for forecast uncertainty in the inventory decisions.

The code demonstrates a production-ready framework that can be adapted to various inventory management scenarios across different industries.