Optimizing Machine Maintenance Planning with Markov Decision Processes (MDP)

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

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

📘 Problem Overview

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

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

Each day, you can choose one of two actions:

  • 0 (Do nothing)
  • 1 (Repair)

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

🧮 MDP Formulation

States:

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

Actions:

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

Transition probabilities:

If you “do nothing”:

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

If you “repair”:

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

Rewards (actually, costs):

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

We’ll minimize these costs using value iteration.

🐍 Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
import numpy as np
import matplotlib.pyplot as plt

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

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

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

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

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

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

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

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

🔍 Code Explanation

Transition Probabilities and Costs

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

Value Iteration

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

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

📊 Visualization

Let’s visualize the convergence and the resulting policy.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# Convergence plot
plt.figure(figsize=(8, 4))
plt.plot(delta_list)
plt.title("Value Iteration Convergence")
plt.xlabel("Iteration")
plt.ylabel("Delta (max change in V)")
plt.grid(True)
plt.show()

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

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

📈 Result Interpretation

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


  • The value iteration converged smoothly, as shown in the delta plot.

  • The optimal value function indicates the minimum expected long-term cost from each state:

    • Good: 48.74
    • Worn: 56.30
    • Broken: 56.30

  • The optimal policy is:

    • In “Good” state: Do nothing
    • In “Worn” state: Repair
    • In “Broken” state: Repair

This policy makes sense:

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

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

✅ Conclusion

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

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

Optimizing Store Layout for Customer Flow and Sales Maximization with Python

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

Let’s dive right in!

Problem Setup

Imagine a small store selling four types of products:

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

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

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

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

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

Mathematical Formulation

Let:

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

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

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

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

Python Code

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np
import itertools
import matplotlib.pyplot as plt

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

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

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

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

# Evaluate expected sales for each arrangement
results = []

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

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

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

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

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

Code Explanation

  • Products and Matrices:
    We defined the four products and created a matrix P representing the probability that a customer moves from one product to another.
    S holds the purchase probability at each product.

  • Arrangements:
    We used itertools.permutations to generate all possible placements of products. (Since there are 4 products, there are $(4! = 24)$ possible layouts.)

  • Sales Evaluation:
    For each arrangement, we calculated the expected sales by summing over all transitions, weighted by transition probability and purchase probability at the destination.

  • Finding the Best Layout:
    We simply picked the arrangement with the maximum expected sales.

  • Visualization:
    Finally, we plotted the expected sales for all possible layouts using a horizontal bar chart. Higher bars indicate better layouts!

Result and Visualization

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

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

Graph:

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

Conclusion

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

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

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

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

Solving the Set Cover Problem with Python

An Optimization Example

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

What is the Set Cover Problem?

The Set Cover Problem can be stated as follows:

Given:

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

Find:

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

In mathematical terms, we want to:

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import numpy as np
import matplotlib.pyplot as plt
import itertools
from matplotlib.patches import Circle
import pulp

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

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

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

# Define coverage radius
coverage_radius = 3.0

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Code Explanation

Let’s break down this implementation step by step:

1. Problem Setup

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

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

2. Computing the Coverage Matrix

We compute a binary coverage matrix where:

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

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

3. Mathematical Formulation

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

Variables:

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

Objective:

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

Constraints:

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

4. Solving with PuLP

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

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

5. Results Visualization

We visualize:

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

6. Solution Analysis

We also analyze:

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

Results

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


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

Results and Insights

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

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

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

Practical Applications

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

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

Conclusion

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

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

Minimizing Staff Requirements in Call Centers with Service Level Constraints

A Practical Python Solution

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

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

The Mathematical Foundation

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

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

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

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

Where:

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

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

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial
import pandas as pd
import seaborn as sns

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Code Explanation

Let me break down the key components of the solution:

  1. Erlang C Formula Implementation:

    • The erlang_c() function calculates the probability that a call will need to wait in a queue
    • It handles the complex math behind queueing theory, including calculating system utilization (rho) and idle probability (P0)

  2. Wait Time Probability:

    • The wait_probability() function extends the Erlang C calculation to determine the probability that a customer waits longer than a specific time threshold
    • This is crucial for evaluating our service level constraint (80% of calls answered within 20 seconds)

  3. Agent Optimization:

    • min_agents_for_service_level() is the core optimization function
    • It starts with the theoretical minimum number of agents (ceil(λ/μ)) needed for a stable system
    • Then it incrementally increases the agent count until the service level requirement is met

  4. Visualization and Analysis:

    • Three different visualizations help understand the problem:
      • Service level vs. number of agents
      • Sensitivity analysis of call volume impact
      • A heatmap showing service level for different agent/call volume combinations

Results Analysis

Minimum agents needed: 13
Service level achieved: 82.70%



For our example call center with:

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

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

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

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

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

Practical Applications

This model can be extremely valuable for:

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

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

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

Solving the Maximum Flow Problem in Python with Visualization

In this post, we’ll explore the maximum flow problem, a classic optimization problem in graph theory.
We’ll walk through a specific example, implement the solution in Python, and visualize the results for better understanding.

🧠 What is the Maximum Flow Problem?

The maximum flow problem asks:

Given a network of nodes connected by edges with certain capacities, what is the greatest possible flow from a source node to a sink node without violating any edge capacities?

We can model many real-world systems using this technique — water pipelines, traffic systems, data flow in networks, etc.

The problem is formally defined on a directed graph $( G = (V, E) )$ with:

  • A source node $( s )$
  • A sink node $( t )$
  • A capacity $( c(u, v) \geq 0 ) for every edge ( (u, v) \in E )$

Our goal is to find a flow $( f(u, v) )$ that satisfies:

  1. Capacity constraint: $( 0 \leq f(u, v) \leq c(u, v) )$
  2. Flow conservation: $( \sum f(u, v) = \sum f(v, w) ) for all ( v \neq s, t )$

We’ll solve this using the Edmonds-Karp algorithm, an implementation of the Ford-Fulkerson method using BFS.

🌉 Sample Network

Let’s consider a network with 6 nodes labeled 0 through 5.

  • Node 0 is the source
  • Node 5 is the sink

Here is our network with capacities:

1
2
3
4
5
6
7
8
9
10
(0) --> (1) capacity 16
(0) --> (2) capacity 13
(1) --> (2) capacity 10
(1) --> (3) capacity 12
(2) --> (1) capacity 4
(2) --> (4) capacity 14
(3) --> (2) capacity 9
(3) --> (5) capacity 20
(4) --> (3) capacity 7
(4) --> (5) capacity 4

🧑‍💻 Python Code (with NetworkX)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import networkx as nx
import matplotlib.pyplot as plt

# Create a directed graph
G = nx.DiGraph()

# Add edges and their capacities
edges = [
    (0, 1, 16),
    (0, 2, 13),
    (1, 2, 10),
    (1, 3, 12),
    (2, 1, 4),
    (2, 4, 14),
    (3, 2, 9),
    (3, 5, 20),
    (4, 3, 7),
    (4, 5, 4)
]

for u, v, c in edges:
    G.add_edge(u, v, capacity=c)

# Compute the maximum flow using Edmonds-Karp algorithm
flow_value, flow_dict = nx.maximum_flow(G, 0, 5, flow_func=nx.algorithms.flow.edmonds_karp)

print("Maximum flow:", flow_value)
print("Flow per edge:")
for u in flow_dict:
    for v in flow_dict[u]:
        if flow_dict[u][v] > 0:
            print(f"  {u} -> {v}: {flow_dict[u][v]}")

🔍 Code Explanation

  • We use networkx.DiGraph() to create a directed graph.
  • Each edge has a capacity, set via the capacity= attribute.
  • nx.maximum_flow() returns two things:
    • flow_value: the total maximum flow from source to sink
    • flow_dict: a dictionary with the flow values per edge

This algorithm internally uses Breadth-First Search (BFS) to find augmenting paths in the residual graph and updates flows iteratively.

📈 Visualizing the Flow

Let’s now visualize the original network and overlay the actual flows after solving the problem.

Maximum flow: 23
Flow per edge:
  0 -> 1: 12
  0 -> 2: 11
  1 -> 3: 12
  2 -> 4: 11
  3 -> 5: 19
  4 -> 3: 7
  4 -> 5: 4


1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Draw the graph with flow values
pos = nx.spring_layout(G, seed=42)  # layout for consistent node positioning
edge_labels = {
    (u, v): f"{flow_dict[u][v]}/{G[u][v]['capacity']}"
    for u, v in G.edges()
}

plt.figure(figsize=(10, 6))
nx.draw_networkx(G, pos, with_labels=True, node_size=800, node_color='lightblue')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw_networkx_edges(G, pos, width=2, edge_color='gray')
plt.title("Maximum Flow in the Network (flow/capacity)")
plt.axis('off')
plt.show()

📊 Interpretation

  • The labels like 10/16 mean 10 units of flow used out of 16 units of capacity on that edge.
  • The final output prints the maximum flow value from node 0 to node 5, which is:

$( \boxed{23} )$

This matches the theoretical solution and demonstrates how efficiently Python and NetworkX can solve and visualize network flow problems.

🚀 Summary

We’ve:

  • Introduced the maximum flow problem
  • Created a network and solved it using Edmonds-Karp algorithm
  • Visualized the results using NetworkX + Matplotlib

This approach is great for teaching, testing, and even small-scale simulations of flow networks.

Inventory Optimization in the Supply Chain

📦 Tackling the Bullwhip Effect with Python

The bullwhip effect is a notorious issue in supply chains, where small fluctuations in customer demand cause progressively larger swings in orders and inventory upstream.
This leads to inefficiencies like overstocking, understocking, and increased costs.

In this article, we’ll explore how to simulate and mitigate this effect through inventory optimization using Python.
We’ll model a simple three-tier supply chain (Retailer → Wholesaler → Manufacturer) and apply a smoothing strategy to reduce the variance of demand signals.

🧠 The Problem: Bullwhip Effect

When demand fluctuates at the retail level, the lack of coordination and communication causes upstream suppliers to overreact. To visualize this:

Let $( D_t )$ be the demand at time $( t )$.
Without a smoothing mechanism, the upstream orders $( O_t )$ are directly based on $( D_t )$, leading to amplified volatility.

We aim to implement a basic inventory policy with exponential smoothing:

Where:

  • $( \hat{D}_t )$ is the forecasted demand,
  • $( \alpha \in [0, 1] )$ is the smoothing parameter.

🛠️ The Simulation Code

Here’s a full Python simulation and visualization in Google Colab:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import numpy as np
import matplotlib.pyplot as plt

# Parameters
np.random.seed(42)
T = 50  # time steps
alpha = 0.3  # smoothing parameter

# Generate customer demand: small random walk
customer_demand = np.round(np.cumsum(np.random.normal(0, 1, T)) + 20).astype(int)
customer_demand = np.clip(customer_demand, 15, 30)

# Initialize arrays
forecast_demand = np.zeros(T)
retailer_orders = np.zeros(T)
wholesaler_orders = np.zeros(T)
manufacturer_orders = np.zeros(T)

# Initial forecast
forecast_demand[0] = customer_demand[0]

# Forecast using exponential smoothing
for t in range(1, T):
    forecast_demand[t] = alpha * customer_demand[t-1] + (1 - alpha) * forecast_demand[t-1]

# Order policies
retailer_orders = forecast_demand + 2  # +2 as safety stock
wholesaler_orders = np.roll(retailer_orders, 1)
manufacturer_orders = np.roll(wholesaler_orders, 1)
wholesaler_orders[0] = retailer_orders[0]
manufacturer_orders[0] = wholesaler_orders[0]

# Plot
plt.figure(figsize=(12, 6))
plt.plot(customer_demand, label='Customer Demand', marker='o')
plt.plot(retailer_orders, label='Retailer Orders', marker='o')
plt.plot(wholesaler_orders, label='Wholesaler Orders', marker='o')
plt.plot(manufacturer_orders, label='Manufacturer Orders', marker='o')
plt.title('Bullwhip Effect Simulation with Smoothing')
plt.xlabel('Time')
plt.ylabel('Order Quantity')
plt.legend()
plt.grid(True)
plt.show()

🧩 Code Breakdown

🔹 Demand Generation

We simulate customer demand using a random walk (with clipping to keep values realistic).

1
2
customer_demand = np.round(np.cumsum(np.random.normal(0, 1, T)) + 20).astype(int)
customer_demand = np.clip(customer_demand, 15, 30)

This simulates a fluctuating but bounded customer demand signal.

🔹 Forecasting with Exponential Smoothing

1
forecast_demand[t] = alpha * customer_demand[t-1] + (1 - alpha) * forecast_demand[t-1]

We use exponential smoothing to estimate the upcoming demand based on recent observations.
The smoothing factor $( \alpha = 0.3 )$ gives moderate weight to recent demand.

🔹 Order Policies

Each player (retailer, wholesaler, manufacturer) bases their orders on forecasts with a simple safety stock adjustment.

1
2
3
retailer_orders = forecast_demand + 2
wholesaler_orders = np.roll(retailer_orders, 1)
manufacturer_orders = np.roll(wholesaler_orders, 1)

The .roll() simulates time delay in receiving downstream orders.

📈 Visualizing the Bullwhip Effect

Here’s what the chart shows:

  • Customer Demand stays relatively smooth.
  • Retailer Orders react slightly, thanks to smoothing.
  • Wholesaler Orders and Manufacturer Orders still show oscillations — but far less severe than without smoothing.

By applying even basic forecasting, the upstream order variability is reduced, and the supply chain becomes more stable.

🧠 Takeaways

  • The bullwhip effect is a real challenge — but simple demand smoothing helps.
  • Forecasting helps reduce the variance and lag in ordering behavior.
  • More advanced methods (like ARIMA, machine learning, or multi-echelon optimization) can further reduce volatility.

Optimizing Pricing for Products and Services

Finding the Sweet Spot

Have you ever wondered how companies determine the perfect price for their products? Too high, and customers flee; too low, and profits vanish.
Today, we’ll explore the fascinating world of price optimization using Python! We’ll build a practical model that helps businesses find that pricing sweet spot across multiple products and services.

Let me walk you through a concrete example of optimizing prices for a tech company offering multiple subscription tiers.
We’ll use advanced optimization techniques, visualize the results, and explain the economics behind the decisions.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
# Importing necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Set the aesthetic style of the plots
plt.style.use('ggplot')
sns.set_palette("colorblind")

# Define the demand function for multiple products
def demand_function(prices, base_demand, price_sensitivity, cross_elasticity_matrix):
    """
    Calculate the demand for each product given their prices and elasticities.
    
    Parameters:
    - prices: Array of prices for each product
    - base_demand: Base demand when price is at reference level
    - price_sensitivity: Own-price elasticity for each product
    - cross_elasticity_matrix: Matrix of cross-price elasticities
    
    Returns:
    - Array of demand values for each product
    """
    # Initialize demand with base demand
    demand = base_demand.copy()
    
    # Apply own-price effects
    for i in range(len(prices)):
        demand[i] *= np.exp(-price_sensitivity[i] * prices[i])
    
    # Apply cross-price effects
    for i in range(len(prices)):
        for j in range(len(prices)):
            if i != j:
                # Positive cross-elasticity means products are substitutes
                # Negative cross-elasticity means products are complements
                demand[i] *= np.exp(cross_elasticity_matrix[i][j] * prices[j])
    
    return demand

# Define the profit function to maximize
def profit_function(prices, base_demand, price_sensitivity, cross_elasticity_matrix, costs):
    """
    Calculate the total profit across all products.
    
    Parameters:
    - prices: Array of prices for each product
    - base_demand: Base demand when price is at reference level
    - price_sensitivity: Own-price elasticity for each product
    - cross_elasticity_matrix: Matrix of cross-price elasticities
    - costs: Variable costs for each product
    
    Returns:
    - Total profit (negative for minimization algorithm)
    """
    demand = demand_function(prices, base_demand, price_sensitivity, cross_elasticity_matrix)
    profit = np.sum((prices - costs) * demand)
    # Return negative profit for minimization (scipy.optimize.minimize minimizes the function)
    return -profit

# Define a function to analyze price elasticity across a range
def price_elasticity_analysis(base_price, product_index, price_range_percentage, 
                            base_demand, price_sensitivity, cross_elasticity_matrix, costs):
    """
    Analyze how profit and demand change as one product's price changes.
    
    Parameters:
    - base_price: Base prices for all products
    - product_index: Index of the product to vary
    - price_range_percentage: Percentage range around base price to analyze
    - Other parameters as defined in previous functions
    
    Returns:
    - DataFrame with price, demand, and profit data
    """
    min_price = base_price[product_index] * (1 - price_range_percentage/100)
    max_price = base_price[product_index] * (1 + price_range_percentage/100)
    price_points = np.linspace(min_price, max_price, 100)
    
    results = []
    for price in price_points:
        # Create a copy of the base price and modify the selected product's price
        current_prices = base_price.copy()
        current_prices[product_index] = price
        
        # Calculate demand and profit
        demand = demand_function(current_prices, base_demand, price_sensitivity, cross_elasticity_matrix)
        profit = -profit_function(current_prices, base_demand, price_sensitivity, 
                                cross_elasticity_matrix, costs)
        
        # Store results
        results.append({
            'Price': price,
            'Demand': demand[product_index],
            'Total Demand': np.sum(demand),
            'Profit': profit
        })
    
    return pd.DataFrame(results)

# Set up an example for a company offering three subscription tiers
# Product 0: Basic tier
# Product 1: Premium tier
# Product 2: Enterprise tier

# Initial parameters
num_products = 3
initial_prices = np.array([10.0, 25.0, 50.0])  # Initial prices
base_demand = np.array([5000, 2000, 1000])     # Base demand at reference prices
costs = np.array([2.0, 5.0, 10.0])             # Variable costs per unit

# Price sensitivity (how demand responds to product's own price)
# Higher values mean demand drops more as price increases
price_sensitivity = np.array([0.15, 0.1, 0.05])

# Cross-price elasticity matrix
# Positive values: substitutes (if price of j increases, demand for i increases)
# Negative values: complements (if price of j increases, demand for i decreases)
cross_elasticity_matrix = np.array([
    [0.0,  0.08, 0.04],  # Impact of products on Basic tier
    [0.03, 0.0,  0.06],  # Impact of products on Premium tier
    [0.01, 0.05, 0.0]    # Impact of products on Enterprise tier
])

# Run the optimization
result = minimize(
    profit_function,
    initial_prices,
    args=(base_demand, price_sensitivity, cross_elasticity_matrix, costs),
    method='SLSQP',
    bounds=[(cost*1.1, cost*10) for cost in costs],  # Price must be at least 10% above cost
    options={'disp': True}
)

# Extract optimized prices
optimal_prices = result.x
print(f"Optimal prices: {optimal_prices}")

# Calculate demand and profit at optimal prices
optimal_demand = demand_function(optimal_prices, base_demand, price_sensitivity, cross_elasticity_matrix)
optimal_profit = -profit_function(optimal_prices, base_demand, price_sensitivity, 
                               cross_elasticity_matrix, costs)

print(f"Demand at optimal prices: {optimal_demand}")
print(f"Total profit at optimal prices: ${optimal_profit:.2f}")

# Analyze how changing each product's price affects overall profit
price_range = 50  # Analyze prices ±50% around optimal
product_names = ["Basic Tier", "Premium Tier", "Enterprise Tier"]

# Create visualizations

# 1. Bar chart comparing initial vs. optimal prices
plt.figure(figsize=(12, 6))
width = 0.35
x = np.arange(num_products)
plt.bar(x - width/2, initial_prices, width, label='Initial Prices')
plt.bar(x + width/2, optimal_prices, width, label='Optimal Prices')
plt.xlabel('Product')
plt.ylabel('Price ($)')
plt.title('Comparison of Initial vs. Optimal Prices')
plt.xticks(x, product_names)
plt.legend()
plt.tight_layout()
plt.savefig('price_comparison.png', dpi=300)
plt.show()

# 2. Plot profit curves for each product around optimal price
plt.figure(figsize=(16, 6))

for i in range(num_products):
    plt.subplot(1, 3, i+1)
    analysis_df = price_elasticity_analysis(
        optimal_prices, i, price_range, base_demand, 
        price_sensitivity, cross_elasticity_matrix, costs
    )
    
    # Find the optimal price within our analyzed range
    max_profit_idx = analysis_df['Profit'].idxmax()
    max_profit_price = analysis_df.loc[max_profit_idx, 'Price']
    
    plt.plot(analysis_df['Price'], analysis_df['Profit'])
    plt.axvline(x=optimal_prices[i], color='red', linestyle='--', 
                label=f'Optimal: ${optimal_prices[i]:.2f}')
    plt.title(f'Profit Curve for {product_names[i]}')
    plt.xlabel('Price ($)')
    plt.ylabel('Total Profit ($)')
    plt.legend()
    plt.grid(True)

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

# 3. Plot demand curves
plt.figure(figsize=(16, 6))

for i in range(num_products):
    plt.subplot(1, 3, i+1)
    analysis_df = price_elasticity_analysis(
        optimal_prices, i, price_range, base_demand, 
        price_sensitivity, cross_elasticity_matrix, costs
    )
    
    plt.plot(analysis_df['Price'], analysis_df['Demand'])
    plt.axvline(x=optimal_prices[i], color='red', linestyle='--', 
                label=f'Optimal: ${optimal_prices[i]:.2f}')
    plt.title(f'Demand Curve for {product_names[i]}')
    plt.xlabel('Price ($)')
    plt.ylabel('Quantity Demanded')
    plt.legend()
    plt.grid(True)

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

# 4. Create a 3D surface plot for two products to visualize the profit landscape
# We'll vary prices for Basic and Premium tiers while keeping Enterprise at optimal
plt.figure(figsize=(12, 10))
ax = plt.axes(projection='3d')

# Create a grid of prices for products 0 and 1
price_range_percent = 30
p0_range = np.linspace(optimal_prices[0] * (1 - price_range_percent/100), 
                     optimal_prices[0] * (1 + price_range_percent/100), 20)
p1_range = np.linspace(optimal_prices[1] * (1 - price_range_percent/100), 
                     optimal_prices[1] * (1 + price_range_percent/100), 20)
P0, P1 = np.meshgrid(p0_range, p1_range)
profit_values = np.zeros(P0.shape)

# Calculate profit for each price combination
for i in range(len(p0_range)):
    for j in range(len(p1_range)):
        current_prices = optimal_prices.copy()
        current_prices[0] = P0[i, j]
        current_prices[1] = P1[i, j]
        profit_values[i, j] = -profit_function(current_prices, base_demand, 
                                             price_sensitivity, cross_elasticity_matrix, costs)

# Create the surface plot
surf = ax.plot_surface(P0, P1, profit_values, cmap=cm.coolwarm,
                       linewidth=0, antialiased=True, alpha=0.8)

# Mark the optimal point
ax.scatter([optimal_prices[0]], [optimal_prices[1]], 
           [-profit_function(optimal_prices, base_demand, price_sensitivity, 
                          cross_elasticity_matrix, costs)],
           color='black', s=100, label='Optimal Price Point')

ax.set_xlabel('Basic Tier Price ($)')
ax.set_ylabel('Premium Tier Price ($)')
ax.set_zlabel('Profit ($)')
ax.set_title('Profit Landscape for Basic and Premium Tiers')
plt.colorbar(surf, ax=ax, shrink=0.5, aspect=5, label='Profit ($)')
plt.savefig('profit_landscape_3d.png', dpi=300)
plt.show()

# 5. Create a heatmap for the same profit landscape
plt.figure(figsize=(10, 8))
profit_df = pd.DataFrame(profit_values, index=p0_range, columns=p1_range)
sns.heatmap(profit_df, cmap='viridis', annot=False)
plt.xlabel('Premium Tier Price ($)')
plt.ylabel('Basic Tier Price ($)')
plt.title('Profit Heatmap for Basic vs Premium Tier Pricing')
plt.tight_layout()
plt.savefig('profit_heatmap.png', dpi=300)
plt.show()

# 6. Create a summary table of results
summary_data = {
    'Product': product_names,
    'Initial Price ($)': initial_prices,
    'Optimal Price ($)': optimal_prices,
    'Price Change (%)': (optimal_prices - initial_prices) / initial_prices * 100,
    'Demand at Optimal': optimal_demand,
    'Unit Profit ($)': optimal_prices - costs,
    'Total Profit ($)': (optimal_prices - costs) * optimal_demand
}

summary_df = pd.DataFrame(summary_data)
summary_df['Price Change (%)'] = summary_df['Price Change (%)'].round(2)
summary_df['Unit Profit ($)'] = summary_df['Unit Profit ($)'].round(2)
summary_df['Total Profit ($)'] = summary_df['Total Profit ($)'].round(2)

print("\nSummary of Pricing Optimization Results:")
print(summary_df)

# Display the cross-elasticity matrix as a heatmap to show relationships between products
plt.figure(figsize=(8, 6))
sns.heatmap(cross_elasticity_matrix, annot=True, cmap='coolwarm', 
            xticklabels=product_names, yticklabels=product_names)
plt.title('Cross-Price Elasticity Matrix')
plt.xlabel('Price Change for Product')
plt.ylabel('Demand Impact on Product')
plt.tight_layout()
plt.savefig('cross_elasticity_heatmap.png', dpi=300)
plt.show()

The Pricing Optimization Model Explained

Our pricing model focuses on finding the optimal prices for three subscription tiers (Basic, Premium, and Enterprise) offered by a tech company.
Let’s break down the key components:

Understanding the Mathematical Framework

The foundation of our pricing model is based on economic principles of demand elasticity.
The demand function is modeled as:

$$D_i(p) = D_i^0 \cdot e^{-\epsilon_i p_i} \cdot \prod_{j \neq i} e^{\epsilon_{ij} p_j}$$

Where:

  • $D_i(p)$ is the demand for product $i$
  • $D_i^0$ is the base demand for product $i$
  • $\epsilon_i$ is the own-price elasticity for product $i$
  • $\epsilon_{ij}$ is the cross-price elasticity between products $i$ and $j$
  • $p_i$ and $p_j$ are the prices of products $i$ and $j$

The profit function we’re maximizing is:

$$\Pi(p) = \sum_{i=1}^{n} (p_i - c_i) \cdot D_i(p)$$

Where $c_i$ is the variable cost for product $i$.

Key Components of the Code

  1. Demand Function: Models how demand changes with price, incorporating both own-price effects (how a product’s price affects its own demand) and cross-price effects (how other products’ prices affect demand).

  2. Profit Function: Calculates total profit across all products based on the demand function, prices, and costs.

  3. Optimization: Uses scipy.optimize.minimize to find the price combination that maximizes profit, with constraints that prices must be at least 10% above cost.

  4. Elasticity Analysis: Examines how changes in each product’s price affect demand and profit.

The Visualization Suite

We’ve created several visualizations to understand the pricing dynamics:

  1. Price Comparison Chart: Shows initial vs. optimal prices for each tier

    Optimization terminated successfully    (Exit mode 0)
         Current function value: -13776186.529477736
         Iterations: 8
         Function evaluations: 18
         Gradient evaluations: 4
Optimal prices: [19.98833796 49.99915445 99.99994941]
Demand at optimal prices: [742205   9278     85]
Total profit at optimal prices: $13776186.53

  2. Profit Curves: Illustrates how profit changes with price for each product

  1. Demand Curves: Shows the relationship between price and quantity demanded

  1. 3D Profit Landscape: Visualizes how profits change when varying two products’ prices simultaneously

  1. Profit Heatmap: A 2D representation of the profit landscape

Summary of Pricing Optimization Results:
           Product  Initial Price ($)  Optimal Price ($)  Price Change (%)  \
0       Basic Tier               10.0          19.988338             99.88   
1     Premium Tier               25.0          49.999154            100.00   
2  Enterprise Tier               50.0          99.999949            100.00   

   Demand at Optimal  Unit Profit ($)  Total Profit ($)  
0             742205            17.99       13351034.38  
1               9278            45.00         417502.16  
2                 85            90.00           7650.00
  1. Cross-Elasticity Heatmap: Shows relationships between products (substitutes vs. complements)

Analysis of Results

When we run the optimization, we find that the optimal prices differ significantly from our initial guesses. Let’s analyze why:

Price-Demand Relationships

The Basic tier shows the highest price sensitivity (0.15), meaning customers are more price-conscious at this level. The Premium tier has moderate sensitivity (0.10), while Enterprise customers are least price-sensitive (0.05), which makes sense as enterprise clients often value features over cost.

Cross-Product Effects

The cross-elasticity matrix reveals interesting product relationships:

  • The Premium tier is a moderate substitute for Basic (0.08)
  • The Enterprise tier strongly influences Premium demand (0.06)
  • Raising the Premium price has a notable effect on Enterprise demand (0.05)

These relationships mean we can’t optimize each product in isolation—we need to consider the entire product ecosystem.

Profit Maximization

The 3D profit landscape reveals that small deviations from optimal pricing can significantly impact profitability. The steep gradient near the peak indicates high sensitivity to price changes, especially for the Basic tier.

The profit heatmap confirms this finding, showing how the combination of Basic and Premium tier pricing creates “sweet spots” of profitability.
Finding this optimal combination is crucial for maximizing overall business value.

Business Implications

The optimal prices we’ve calculated have several important business implications:

  1. Price Differentiation: The significant price gaps between tiers help segment the market effectively.

  2. Cross-Selling Opportunities: Understanding cross-elasticities reveals opportunities for targeted promotions and bundles.

  3. Margin Management: The unit profit for each tier shows where the business makes most of its money, which can guide feature development priorities.

  4. Competitive Positioning: The model can be extended to incorporate competitor prices and market share dynamics.

Extending the Model

This pricing model can be extended in several ways:

  1. Time-Series Analysis: Incorporate seasonal demand patterns
  2. Customer Segmentation: Different elasticities for different customer groups
  3. Feature-Based Pricing: Break down subscription tiers into feature components
  4. Dynamic Pricing: Update prices based on changing market conditions

Conclusion

Optimal pricing is both art and science.
Our model provides a data-driven foundation for pricing decisions while giving business leaders flexibility to incorporate qualitative factors.
By understanding the mathematical relationships between prices, demand, and profits across multiple products, companies can make better pricing decisions that boost both revenues and customer satisfaction.

Optimizing Production Planning Based on Demand Forecasting Using Python

In today’s data-driven manufacturing environments, aligning production with future demand is crucial.
Producing too much leads to high inventory costs; producing too little results in lost sales.
In this post, we’ll use Python to forecast demand and build an optimal production plan that minimizes cost while satisfying predicted demand.

Let’s dive into a real-world-inspired example!

🧠 Problem Overview

Imagine a factory producing a single product.
Based on historical sales data, we forecast future monthly demand.
The goal is to decide how much to produce each month to minimize total cost.

Given:

  • Forecasted demand for the next 6 months
  • Production cost per unit: $10
  • Inventory holding cost per unit per month: $2
  • Initial inventory: 0
  • Maximum monthly production capacity: 100 units

Objective:

Minimize the total cost = production cost + inventory cost
Subject to meeting the demand in each month.

Let:

  • $( d_t )$: forecasted demand in month $( t )$
  • $( x_t )$: units produced in month $( t )$
  • $( I_t )$: inventory at the end of month $( t )$

The constraints:

  • Inventory balance:
    $$
    I_t = I_{t-1} + x_t - d_t
    $$
  • Non-negativity:
    $$
    x_t \geq 0,\quad I_t \geq 0
    $$
  • Production capacity:
    $$
    x_t \leq 100
    $$

🧮 Forecasted Demand

Let’s assume we have this demand forecast for the next 6 months:

1
2
3
4
5
import numpy as np
import matplotlib.pyplot as plt

months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun']
demand = np.array([80, 100, 60, 90, 110, 70])

🧑‍💻 Optimization with PuLP

We’ll use the PuLP library to solve this as a linear programming problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
!pip install pulp

from pulp import LpMinimize, LpProblem, LpVariable, lpSum, LpStatus

# Problem definition
model = LpProblem("Production_Planning", LpMinimize)

n_months = len(demand)
production = [LpVariable(f"x_{t}", lowBound=0, upBound=100) for t in range(n_months)]
inventory = [LpVariable(f"I_{t}", lowBound=0) for t in range(n_months)]

# Objective function: production cost + inventory cost
production_cost = lpSum(10 * production[t] for t in range(n_months))
inventory_cost = lpSum(2 * inventory[t] for t in range(n_months))
model += production_cost + inventory_cost

# Inventory balance constraints
for t in range(n_months):
    if t == 0:
        model += inventory[t] == production[t] - demand[t]
    else:
        model += inventory[t] == inventory[t-1] + production[t] - demand[t]

# Solve
model.solve()

# Output
print(f"Status: {LpStatus[model.status]}")
for t in range(n_months):
    print(f"{months[t]} - Production: {production[t].value():.0f}, Inventory: {inventory[t].value():.0f}")
Status: Optimal
Jan - Production: 80, Inventory: 0
Feb - Production: 100, Inventory: 0
Mar - Production: 60, Inventory: 0
Apr - Production: 100, Inventory: 10
May - Production: 100, Inventory: 0
Jun - Production: 70, Inventory: 0

📊 Visualization of Results

Now let’s visualize the optimal plan:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
prod_vals = [production[t].value() for t in range(n_months)]
inv_vals = [inventory[t].value() for t in range(n_months)]

plt.figure(figsize=(10, 5))
plt.plot(months, demand, label="Demand", marker='o')
plt.plot(months, prod_vals, label="Production", marker='s')
plt.plot(months, inv_vals, label="Inventory", marker='^')
plt.title("Optimal Production Plan")
plt.xlabel("Month")
plt.ylabel("Units")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

📈 Detailed Explanation

🛠 Objective Function

The cost function:
$$
\text{Total Cost} = \sum_{t=1}^{6} \left(10 \cdot x_t + 2 \cdot I_t \right)
$$
It combines the per-unit production and inventory costs.

📦 Inventory Balance

Each month’s ending inventory is the previous month’s inventory plus new production minus demand:
$$
I_t = I_{t-1} + x_t - d_t
$$

🔍 Results Interpretation

From the printed results and graph:

  • Production is smoothed to avoid overproduction.
  • Inventory is used strategically to satisfy high demand in future months.
  • Costs are minimized by balancing production with holding costs.

🧾 Conclusion

This example illustrates how even a basic demand forecast can drive powerful optimizations in production planning.

By using Python and PuLP, businesses can automate and optimize their manufacturing strategies, saving costs and improving efficiency.

Zone-Based Delivery Optimization with Clustering and TSP in Python

🚚 Optimizing Zone-Based Delivery: Clustering Meets Route Optimization in Python

When planning delivery operations across a city, efficiency is everything.
One common approach is zone delivery planning, where delivery addresses are grouped into zones, and then an optimal route is calculated within each zone.
This strategy reduces travel time, improves scheduling, and enhances customer satisfaction.

In this post, we’ll walk through a practical example using Python where we:

  1. Generate delivery points (randomly for simulation).
  2. Cluster them into zones using K-Means clustering.
  3. Use Google OR-Tools to find the optimal delivery route inside each zone.
  4. Visualize everything step-by-step.

Let’s get into it!

🧠 Step 1: Simulating Delivery Points

We start by creating 50 random delivery points in a city-like area.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np
import matplotlib.pyplot as plt

# Generate random delivery locations
np.random.seed(42)
num_points = 50
X = np.random.rand(num_points, 2) * 100  # Coordinates within 100x100 grid

# Plot delivery points
plt.figure(figsize=(6, 6))
plt.scatter(X[:, 0], X[:, 1], c='blue')
plt.title("Delivery Points")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
plt.show()

📦 Step 2: Clustering Into Zones with K-Means

Now, we’ll use K-Means to split these delivery points into, say, 4 zones.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from sklearn.cluster import KMeans

num_clusters = 4
kmeans = KMeans(n_clusters=num_clusters, random_state=42)
labels = kmeans.fit_predict(X)

# Plot clustered points
plt.figure(figsize=(6, 6))
for i in range(num_clusters):
    cluster = X[labels == i]
    plt.scatter(cluster[:, 0], cluster[:, 1], label=f'Zone {i+1}')
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], 
            c='black', marker='x', s=100, label='Centers')
plt.title("Clustered Delivery Zones")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.legend()
plt.grid(True)
plt.show()

🗺️ Step 3: Route Optimization with Google OR-Tools

For each zone, we’ll use OR-Tools to solve the Traveling Salesman Problem (TSP).

TSP Mathematical Formulation

Given a set of nodes $ ( V = {v_1, v_2, …, v_n} )$, and distances $( d_{ij} )$, the goal is to find the shortest possible route that visits each node once and returns to the starting point:

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

Subject to:

$$
\sum_{j=1}^{n} x_{ij} = 1 \quad \forall i,\quad
\sum_{i=1}^{n} x_{ij} = 1 \quad \forall j
$$

Let’s implement this using OR-Tools.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from ortools.constraint_solver import pywrapcp, routing_enums_pb2

def create_distance_matrix(locations):
    from scipy.spatial.distance import cdist
    return cdist(locations, locations).astype(int)

def solve_tsp(distance_matrix):
    manager = pywrapcp.RoutingIndexManager(len(distance_matrix), 1, 0)
    routing = pywrapcp.RoutingModel(manager)

    def distance_callback(from_idx, to_idx):
        return distance_matrix[manager.IndexToNode(from_idx)][manager.IndexToNode(to_idx)]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    search_params = pywrapcp.DefaultRoutingSearchParameters()
    search_params.first_solution_strategy = routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC

    solution = routing.SolveWithParameters(search_params)
    
    if not solution:
        return None
    
    # Extract the route
    index = routing.Start(0)
    route = []
    while not routing.IsEnd(index):
        route.append(manager.IndexToNode(index))
        index = solution.Value(routing.NextVar(index))
    route.append(manager.IndexToNode(index))
    return route

🔁 Solving for Each Zone and Visualizing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
zone_routes = []
plt.figure(figsize=(8, 8))

colors = ['r', 'g', 'b', 'orange']
for zone_id in range(num_clusters):
    zone_points = X[labels == zone_id]
    if len(zone_points) < 2:
        continue  # Skip small clusters

    distance_matrix = create_distance_matrix(zone_points)
    route = solve_tsp(distance_matrix)
    if route is None:
        print(f"No route found for zone {zone_id}")
        continue

    # Save and plot
    zone_routes.append(route)
    zone_route_coords = zone_points[route]

    plt.plot(zone_route_coords[:, 0], zone_route_coords[:, 1], 
             marker='o', color=colors[zone_id], label=f'Zone {zone_id+1}')

plt.title("Optimized Delivery Routes per Zone")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
plt.legend()
plt.show()

📊 Analysis

  • Clustering reduces complexity: Instead of solving TSP for 50 points (computationally expensive), we solve 4 smaller TSPs.
  • Flexibility: Number of clusters can be adapted depending on truck capacity, time windows, or depot proximity.
  • Scalability: This method works for hundreds of points when clusters are balanced.

✅ Conclusion

Zone-based delivery planning is a powerful approach to optimizing logistics.
By combining unsupervised learning (clustering) with combinatorial optimization (TSP), we achieve efficient and scalable routing.

A Python Adventure in Project Management

Finding the Critical Path

Today, I’m excited to share a practical guide to solving Critical Path Method (CPM) problems with Python.
CPM is a fundamental project management technique that helps identify the sequence of activities that determine the shortest time to complete a project.

Let’s dive into a concrete example and see how we can implement the CPM algorithm in Python!

The Problem: Building a New Mobile App

Imagine you’re managing the development of a new mobile app.
The project consists of several activities with dependencies between them.
Your goal is to determine the earliest possible completion time and identify the critical path - the sequence of activities that must not be delayed to finish the project on time.

Here’s our project breakdown:

Activity Description Duration (days) Predecessors
A Requirements gathering 3 None
B System design 4 A
C Database setup 2 A
D Frontend development 5 B
E Backend development 6 B, C
F Integration 3 D, E
G Testing 4 F
H Deployment 2 G

Now, let’s solve this using Python!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Rectangle
from matplotlib.colors import LinearSegmentedColormap

# Define the project activities with their durations and dependencies
activities = {
    'A': {'name': 'Requirements gathering', 'duration': 3, 'predecessors': []},
    'B': {'name': 'System design', 'duration': 4, 'predecessors': ['A']},
    'C': {'name': 'Database setup', 'duration': 2, 'predecessors': ['A']},
    'D': {'name': 'Frontend development', 'duration': 5, 'predecessors': ['B']},
    'E': {'name': 'Backend development', 'duration': 6, 'predecessors': ['B', 'C']},
    'F': {'name': 'Integration', 'duration': 3, 'predecessors': ['D', 'E']},
    'G': {'name': 'Testing', 'duration': 4, 'predecessors': ['F']},
    'H': {'name': 'Deployment', 'duration': 2, 'predecessors': ['G']}
}

# Create a directed graph
G = nx.DiGraph()

# Add nodes (activities) to the graph
for activity_id, activity_info in activities.items():
    G.add_node(activity_id, **activity_info)

# Add edges (dependencies) to the graph
for activity_id, activity_info in activities.items():
    for predecessor in activity_info['predecessors']:
        G.add_edge(predecessor, activity_id)

# Calculate earliest start and finish times (forward pass)
earliest_start = {'A': 0}  # Activity A starts at time 0
earliest_finish = {}

# Topological sort ensures we process activities in the correct order
for activity_id in nx.topological_sort(G):
    activity_info = G.nodes[activity_id]
    
    # Set earliest start time
    if activity_id not in earliest_start:
        if not activity_info['predecessors']:
            earliest_start[activity_id] = 0
        else:
            earliest_start[activity_id] = max(earliest_finish.get(pred, 0) for pred in activity_info['predecessors'])
    
    # Calculate earliest finish time
    earliest_finish[activity_id] = earliest_start[activity_id] + activity_info['duration']

# Calculate latest start and finish times (backward pass)
project_duration = max(earliest_finish.values())
latest_finish = {activity_id: project_duration for activity_id in activities}
latest_start = {}

# Process activities in reverse topological order
for activity_id in reversed(list(nx.topological_sort(G))):
    # Find successors of the current activity
    successors = list(G.successors(activity_id))
    
    if successors:
        latest_finish[activity_id] = min(latest_start.get(succ, project_duration) for succ in successors)
    
    latest_start[activity_id] = latest_finish[activity_id] - G.nodes[activity_id]['duration']

# Calculate float (slack) for each activity
float_times = {}
for activity_id in activities:
    float_times[activity_id] = latest_start[activity_id] - earliest_start[activity_id]

# Identify the critical path (activities with zero float)
critical_path = [activity_id for activity_id, float_time in float_times.items() if float_time == 0]

# Print the results
print(f"Project Duration: {project_duration} days")
print(f"Critical Path: {' -> '.join(critical_path)}")
print("\nActivity Details:")
print(f"{'Activity':<10}{'Description':<25}{'Duration':<10}{'ES':<5}{'EF':<5}{'LS':<5}{'LF':<5}{'Float':<5}{'Critical':<8}")
print("-" * 80)

for activity_id, activity_info in sorted(activities.items()):
    is_critical = activity_id in critical_path
    print(f"{activity_id:<10}{activity_info['name']:<25}{activity_info['duration']:<10}{earliest_start[activity_id]:<5}"
          f"{earliest_finish[activity_id]:<5}{latest_start[activity_id]:<5}{latest_finish[activity_id]:<5}"
          f"{float_times[activity_id]:<5}{'Yes' if is_critical else 'No':<8}")

# Visualization of the network diagram
plt.figure(figsize=(14, 8))
pos = nx.spring_layout(G, seed=42)  # Layout for the graph

# Draw the network
node_colors = ['red' if node in critical_path else 'skyblue' for node in G.nodes()]
edge_colors = ['red' if u in critical_path and v in critical_path else 'black' for u, v in G.edges()]
edge_width = [2 if u in critical_path and v in critical_path else 1 for u, v in G.edges()]

nx.draw(G, pos, with_labels=True, node_color=node_colors, edge_color=edge_colors, 
        width=edge_width, node_size=700, font_size=10, font_weight='bold')

plt.title('Project Network Diagram with Critical Path (in red)', fontsize=16)
plt.savefig('network_diagram.png', dpi=300, bbox_inches='tight')
plt.close()

# Create a Gantt chart
plt.figure(figsize=(14, 8))

# Create custom colormap: blue for normal activities, red for critical path
cmap = LinearSegmentedColormap.from_list('custom_cmap', ['skyblue', 'red'])

# Sort activities by earliest start time
sorted_activities = sorted(activities.items(), key=lambda x: earliest_start[x[0]])

# Draw the Gantt chart
y_positions = np.arange(len(activities))
y_labels = [f"{activity_id}: {info['name']}" for activity_id, info in sorted_activities]

# Draw activity bars
for i, (activity_id, info) in enumerate(sorted_activities):
    is_critical = activity_id in critical_path
    color = 'red' if is_critical else 'skyblue'
    
    # Draw the main activity bar
    plt.barh(i, info['duration'], left=earliest_start[activity_id], color=color, 
             edgecolor='black', alpha=0.8)
    
    # Draw the float/slack time (if any)
    if float_times[activity_id] > 0:
        plt.barh(i, float_times[activity_id], left=earliest_finish[activity_id], 
                 color='lightgray', alpha=0.5, edgecolor='gray', hatch='/')
    
    # Add text labels on the bars
    plt.text(earliest_start[activity_id] + info['duration']/2, i, 
             f"{activity_id} ({info['duration']}d)", 
             ha='center', va='center', color='black', fontweight='bold')

# Draw time grid
for t in range(0, project_duration + 1):
    plt.axvline(x=t, color='gray', linestyle='--', alpha=0.3)

# Set chart properties
plt.yticks(y_positions, y_labels)
plt.xlabel('Time (days)')
plt.grid(axis='x', alpha=0.3)
plt.title('Project Gantt Chart with Critical Path', fontsize=16)

# Add legend
legend_elements = [
    Rectangle((0, 0), 1, 1, color='red', alpha=0.8, label='Critical Activity'),
    Rectangle((0, 0), 1, 1, color='skyblue', alpha=0.8, label='Non-Critical Activity'),
    Rectangle((0, 0), 1, 1, color='lightgray', alpha=0.5, hatch='/', label='Float Time')
]
plt.legend(handles=legend_elements, loc='upper right')

# Show the Gantt chart
plt.tight_layout()
plt.savefig('gantt_chart.png', dpi=300, bbox_inches='tight')
plt.show()

Understanding the Code: Step by Step

Let’s break down the CPM implementation:

1. Setting Up the Project Model

First, we define our project activities with their durations and dependencies.
We use a dictionary structure where each key is an activity ID, and the value contains relevant information about that activity.

We then create a directed graph using NetworkX, a powerful Python library for graph theory.
Each node represents an activity, and edges represent dependencies between activities.

2. Forward Pass: Calculating Earliest Times

The forward pass calculates the earliest possible start and finish times for each activity:

  • Earliest Start Time (ES): The earliest time an activity can begin, which is the maximum of the earliest finish times of all its predecessors.
  • Earliest Finish Time (EF): ES + activity duration.

We use a topological sort to ensure we process activities in the correct order (predecessors before successors).

3. Backward Pass: Calculating Latest Times

The backward pass calculates the latest allowable start and finish times for each activity:

  • Latest Finish Time (LF): The latest time an activity can end without delaying the project, which is the minimum of the latest start times of all its successors.
  • Latest Start Time (LS): LF - activity duration.

4. Calculating Float and Identifying the Critical Path

The float (or slack) for each activity is calculated as LS - ES.
Activities with zero float are on the critical path - these activities must be completed on time to avoid delaying the entire project.

5. Visualization

We create two visualizations:

  1. A network diagram showing the dependencies between activities, with the critical path highlighted in red.
  2. A Gantt chart showing the timeline of activities, distinguishing between critical and non-critical activities, and showing float times.

Mathematical Representation

In CPM, we can represent the timing calculations using the following equations:

For an activity $i$:

  • Earliest Start Time: $ES_i = \max{EF_j | j \in \text{predecessors of } i}$
  • Earliest Finish Time: $EF_i = ES_i + D_i$
  • Latest Finish Time: $LF_i = \min{LS_j | j \in \text{successors of } i}$
  • Latest Start Time: $LS_i = LF_i - D_i$
  • Float Time: $Float_i = LS_i - ES_i$

Where $D_i$ is the duration of activity $i$.

Results and Analysis

When we run our code, we get the following output:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Project Duration: 22 days
Critical Path: A -> B -> E -> F -> G -> H

Activity Details:
Activity  Description                Duration  ES   EF   LS   LF   Float Critical
--------------------------------------------------------------------------------
A         Requirements gathering     3         0    3    0    3    0    Yes     
B         System design              4         3    7    3    7    0    Yes     
C         Database setup             2         3    5    5    7    2    No      
D         Frontend development       5         7    12   10   15   3    No      
E         Backend development        6         7    13   7    13   0    Yes     
F         Integration                3         13   16   13   16   0    Yes     
G         Testing                    4         16   20   16   20   0    Yes     
H         Deployment                 2         20   22   20   22   0    Yes

This tells us:

  1. The project will take 22 days to complete.
  2. The critical path is A → B → E → F → G → H.
  3. Activities C and D have float times of 2 and 3 days respectively, meaning they can be delayed by that much without affecting the overall project duration.

Visualizing the Results

Our code generates two helpful visualizations:

1. Network Diagram

The network diagram shows the relationships between activities.
The critical path is highlighted in red, making it easy to identify the activities that require special attention.

2. Gantt Chart

The Gantt chart provides a timeline view of the project. It shows:

  • When each activity is scheduled to start and finish
  • Which activities are on the critical path (in red)
  • Float time for non-critical activities (gray hatched areas)

This visualization is particularly useful for project managers to track progress and identify potential bottlenecks.

Key Insights from CPM Analysis

  1. Critical Activities: The activities on the critical path (A, B, E, F, G, H) require careful management, as any delay will extend the project completion time.

  2. Resource Allocation: Non-critical activities (C, D) have some flexibility in their scheduling, which allows for better resource allocation.

  3. Potential Optimizations: If we want to shorten the project duration, we should focus on reducing the duration of activities on the critical path.

Conclusion

The Critical Path Method is a powerful tool for project management that helps identify the sequence of activities that determine the minimum time needed to complete a project.
By implementing CPM in Python, we can quickly analyze complex projects and make informed decisions about scheduling and resource allocation.

The next time you’re planning a project, remember to identify your critical path - it might be the key to delivering on time!