Optimizing Power Supply and Demand Balance in Smart Grid Operations

Posted on April 17, 2025 by SmartGridEnthusiast

Smart grids represent the future of electricity distribution, offering intelligent ways to balance power supply and demand.
Today, I’ll walk you through a practical example of how we can optimize this balance using Python and mathematical programming.

Let’s tackle a realistic scenario: managing multiple power sources (traditional and renewable) while meeting consumer demand throughout the day, all while minimizing costs and ensuring stability.

The Smart Grid Optimization Problem

In our example, we’ll model a smart grid with:

  • Multiple power generators (coal, gas, solar, wind)
  • Varying demand throughout a 24-hour period
  • Different production costs and constraints
  • Energy storage capabilities

The mathematical formulation of our problem involves minimizing the total cost while satisfying all constraints:

$$ \min \sum_{t=1}^{T} \sum_{g=1}^{G} C_g \cdot P_{g,t} + \sum_{t=1}^{T} C_{storage} \cdot |S_t - S_{t-1}| $$

Subject to:
$$ \sum_{g=1}^{G} P_{g,t} + S_{t-1} - S_t = D_t \quad \forall t \in T $$
$$ P_{g,t}^{min} \leq P_{g,t} \leq P_{g,t}^{max} \quad \forall g \in G, t \in T $$
$$ 0 \leq S_t \leq S^{max} \quad \forall t \in T $$

Where:

  • $P_{g,t}$ is the power output of generator $g$ at time $t$
  • $C_g$ is the cost coefficient for generator $g$
  • $S_t$ is the stored energy at time $t$
  • $D_t$ is the demand at time $t$

Let’s implement this model using Python with the PuLP library for optimization:

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
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pulp
from matplotlib.ticker import MaxNLocator
import seaborn as sns

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

# Define the problem parameters
num_hours = 24
num_generators = 4

# Generator types and their properties
generator_types = ['Coal', 'Natural Gas', 'Solar', 'Wind']
generator_costs = {'Coal': 100, 'Natural Gas': 80, 'Solar': 5, 'Wind': 10}  # Cost per MWh
generator_min_output = {'Coal': 50, 'Natural Gas': 20, 'Solar': 0, 'Wind': 0}  # Min output in MW
generator_max_output = {'Coal': 200, 'Natural Gas': 150, 'Solar': 100, 'Wind': 120}  # Max output in MW

# Set up time-dependent availability for renewables (solar and wind)
solar_availability = np.zeros(num_hours)
wind_availability = np.zeros(num_hours)

# Solar availability follows a bell curve (daylight hours)
for t in range(num_hours):
    if 6 <= t < 18:  # Daylight hours (6 AM to 6 PM)
        solar_availability[t] = 0.8 * np.sin(np.pi * (t - 6) / 12) + 0.1
    else:
        solar_availability[t] = 0.1  # Minimal solar at night
        
# Wind availability has some variability throughout the day
wind_availability = 0.4 + 0.2 * np.sin(np.pi * np.arange(num_hours) / 12) + 0.1 * np.random.rand(num_hours)

# Demand profile throughout the day (in MW)
demand = 200 + 100 * np.sin(np.pi * (np.arange(num_hours) - 4) / 12) + 50 * np.sin(np.pi * (np.arange(num_hours) - 8) / 8)
demand = demand.astype(int)

# Storage parameters
storage_capacity = 150  # Maximum storage capacity in MWh
storage_cost = 15  # Cost per MWh of charging/discharging
initial_storage = 50  # Initial storage level in MWh

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

# Decision variables
# Power output for each generator at each hour
power_output = {(g, t): pulp.LpVariable(f"Power_{g}_{t}", 
                                         lowBound=generator_min_output[g], 
                                         upBound=generator_max_output[g] * (1 if g not in ['Solar', 'Wind'] else 
                                                                           (solar_availability[t] if g == 'Solar' else wind_availability[t])), 
                                         cat='Continuous')
                for g in generator_types for t in range(num_hours)}

# Storage level at each hour
storage = {t: pulp.LpVariable(f"Storage_{t}", lowBound=0, upBound=storage_capacity, cat='Continuous') 
          for t in range(num_hours)}

# Variables to track charging and discharging (for cost calculation)
charging = {t: pulp.LpVariable(f"Charging_{t}", lowBound=0, cat='Continuous') 
          for t in range(num_hours)}
discharging = {t: pulp.LpVariable(f"Discharging_{t}", lowBound=0, cat='Continuous') 
          for t in range(num_hours)}

# Objective function: minimize total cost
model += pulp.lpSum([generator_costs[g] * power_output[g, t] for g in generator_types for t in range(num_hours)]) + \
         pulp.lpSum([storage_cost * (charging[t] + discharging[t]) for t in range(num_hours)])

# Constraints
# Initial storage level
model += storage[0] == initial_storage

# Storage dynamics
for t in range(1, num_hours):
    model += storage[t] - storage[t-1] == charging[t] - discharging[t]
    
# Power balance: supply = demand at each hour
for t in range(num_hours):
    model += pulp.lpSum([power_output[g, t] for g in generator_types]) + discharging[t] - charging[t] == demand[t]

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

print(f"Status: {pulp.LpStatus[model.status]}")
print(f"Total Cost: ${pulp.value(model.objective):.2f}")

# Extract results
results = pd.DataFrame(index=range(num_hours), columns=generator_types + ['Storage', 'Demand'])

for t in range(num_hours):
    for g in generator_types:
        results.loc[t, g] = power_output[g, t].value()
    results.loc[t, 'Storage'] = storage[t].value()
    results.loc[t, 'Demand'] = demand[t]

# Plot the results
plt.figure(figsize=(15, 10))

# Plot 1: Generator outputs and demand
ax1 = plt.subplot(2, 1, 1)
bottom = np.zeros(num_hours)

for gen in generator_types:
    plt.bar(range(num_hours), results[gen], bottom=bottom, label=gen, alpha=0.7)
    bottom += results[gen].values

plt.plot(range(num_hours), results['Demand'], 'r-', linewidth=2, label='Demand')
plt.title('Hourly Power Generation by Source vs Demand', fontsize=14)
plt.xlabel('Hour of Day')
plt.ylabel('Power (MW)')
plt.legend()
plt.xticks(range(0, num_hours, 2))
plt.grid(True, alpha=0.3)

# Plot 2: Storage level
ax2 = plt.subplot(2, 1, 2)
plt.plot(range(num_hours), results['Storage'], 'g-', linewidth=2, marker='o')
plt.fill_between(range(num_hours), results['Storage'], alpha=0.3, color='green')
plt.title('Energy Storage Level Over Time', fontsize=14)
plt.xlabel('Hour of Day')
plt.ylabel('Storage Level (MWh)')
plt.xticks(range(0, num_hours, 2))
plt.grid(True, alpha=0.3)
plt.tight_layout()

# Additional plot: Cost breakdown
costs = pd.DataFrame(index=generator_types + ['Storage'])
for g in generator_types:
    costs.loc[g, 'Cost'] = sum(generator_costs[g] * power_output[g, t].value() for t in range(num_hours))
    
storage_total_cost = 0
for t in range(num_hours):
    if t > 0:
        storage_total_cost += storage_cost * abs(storage[t].value() - storage[t-1].value())
        
costs.loc['Storage', 'Cost'] = storage_total_cost

plt.figure(figsize=(10, 6))
sns.barplot(x=costs.index, y='Cost', data=costs)
plt.title('Total Cost by Energy Source', fontsize=14)
plt.xlabel('Source')
plt.ylabel('Cost ($)')
plt.grid(True, alpha=0.3)
plt.tight_layout()

# Export results to CSV
results.to_csv('smart_grid_optimization_results.csv')

plt.show()

# Print hourly breakdown
print("\nHourly Power Generation and Storage:")
print(results.round(2))

# Calculate some key metrics
total_renewable = results['Solar'].sum() + results['Wind'].sum()
total_conventional = results['Coal'].sum() + results['Natural Gas'].sum()
total_generation = total_renewable + total_conventional

print(f"\nPercentage of Renewable Energy: {100 * total_renewable / total_generation:.2f}%")
print(f"Maximum Storage Used: {results['Storage'].max():.2f} MWh")
print(f"Storage Utilization: {100 * results['Storage'].max() / storage_capacity:.2f}%")

Code Explanation

Let’s break down this implementation step by step:

1. Problem Setup

First, we define our smart grid parameters:

  • Four generator types (Coal, Natural Gas, Solar, Wind)
  • Cost structures for each generator
  • Operational constraints (min/max output)
  • Time-dependent availability for renewables
  • Demand curve throughout the day
  • Energy storage capabilities

The solar availability follows a realistic daily pattern with peak availability during mid-day hours, while wind power has some natural variability.
Demand follows a common pattern with morning and evening peaks.

2. Mathematical Model Implementation Using PuLP

We use the PuLP library to define our linear programming model:

  • Decision Variables:

    • power_output[g,t]: How much power each generator produces each hour
    • storage[t]: Battery storage level at each hour
    • charging[t] and discharging[t]: Energy flow in/out of storage

  • Objective Function:
    Minimize the total cost of generation plus storage operations.

  • Constraints:

    • Power balance at each hour (supply = demand)
    • Generator output limitations
    • Storage capacity constraints
    • Storage dynamics (tracking energy flow)

3. Solving the Model

The PuLP solver (CBC) finds the optimal solution that minimizes costs while meeting all constraints.
This shows how different power sources should be dispatched throughout the day.

4. Visualization and Analysis

We generate comprehensive visualizations to understand the solution:

  • A stacked bar chart showing how each generator contributes hourly
  • A line plot tracking storage levels
  • A cost breakdown by energy source

Results Analysis

Let’s examine the graphs to understand the optimal solution:

Status: Optimal
Total Cost: $251791.35



Hourly Power Generation and Storage:
    Coal  Natural Gas  Solar   Wind  Storage  Demand
0   50.0        20.00  10.00  33.00    50.00     113
1   50.0        20.00  10.00  65.62    85.62     110
2   50.0        20.00  10.00  68.78   120.40     114
3   50.0        20.00  10.00  72.15   145.56     127
4   50.0        20.00  10.00  70.66   146.22     150
5   50.0        20.00  10.00  73.05   120.27     179
6   50.0        20.00  10.00  72.70    58.97     214
7   50.0        29.75  30.71  81.58     0.00     251
8   50.0       110.00  50.00  76.00     0.00     286
9   50.0       124.96  66.57  73.47     0.00     315
10  50.0       145.47  79.28  60.25     0.00     335
11  50.0       138.88  87.27  65.85     0.00     342
12  50.0       138.01  90.00  57.99     0.00     336
13  50.0       134.39  87.27  44.34     0.00     316
14  50.0       117.54  79.28  38.18     0.00     285
15  50.0        95.20  66.57  33.23     0.00     245
16  50.0        69.13  50.00  30.87     0.00     200
17  50.0        42.18  30.71  31.11     0.00     154
18  50.0        24.82  10.00  29.18     0.00     114
19  50.0        20.00  10.00   3.00     0.00      83
20  50.0        20.00   0.00   0.00     7.00      63
21  50.0        20.00   0.00   0.00    20.00      57
22  50.0        20.00   0.00   0.00    26.00      64
23  50.0        20.00  10.00   4.00    26.00      84

Percentage of Renewable Energy: 42.16%
Maximum Storage Used: 146.22 MWh
Storage Utilization: 97.48%

Generator Dispatch Pattern

Looking at the hourly power generation chart, we can observe several key insights:

  1. Base Load Operation: Coal plants, with their higher minimum output requirements but high operational costs, are used primarily as base load during peak demand hours.

  2. Mid-merit Operation: Natural gas generators, with more flexibility and moderate costs, are ramped up and down to follow the demand curve.

  3. Renewable Integration: Solar energy is maximized during daylight hours, while wind power is used whenever available, as they have the lowest operational costs.

  4. Peak Management: The combination of all sources plus storage is optimized to meet peak demand periods in the morning and evening.

Storage Utilization

The storage level chart shows:

  1. Charge/Discharge Cycles: The storage system charges during low-demand periods (especially when renewable generation is high) and discharges during peak demand periods.

  2. Peak Shaving: The storage helps “shave” demand peaks by providing additional power, reducing the need for expensive peaking generators.

  3. Renewable Smoothing: Storage helps balance the intermittent nature of renewable sources, storing excess generation for later use.

Cost Distribution

The cost breakdown reveals:

  1. Generation Mix Economics: Despite their lower capacity factors, renewables contribute significantly to cost savings.

  2. Storage Economics: The cost of storage operations is offset by the system’s ability to shift energy from low-cost generation periods to high-demand periods.

Practical Applications

This optimization approach can be applied in several real-world scenarios:

  1. Microgrid Management: Isolated communities can optimize their local generation mix.

  2. Utility-Scale Operations: Power companies can optimize their generation fleet.

  3. Renewable Integration Planning: Grid operators can determine optimal storage capacity to accommodate increasing renewable penetration.

  4. Demand Response Programs: Utilities can quantify the value of flexible demand.

Conclusion

Our Python model demonstrates how mathematical optimization can effectively balance power supply and demand in a smart grid environment.
By accounting for the specific characteristics of different generation sources, demand patterns, and storage capabilities, grid operators can minimize costs while ensuring reliable electricity supply.

As renewable energy sources continue to grow in importance, these optimization techniques will become increasingly valuable for managing the complex dynamics of modern power systems.

Real Options Theory in Investment Decision Making

A Practical Example

Today I’m going to walk through a fascinating application of real options theory to investment decision making.
Unlike traditional NPV analysis, real options theory recognizes the value of flexibility in business decisions - the ability to delay, expand, contract or abandon projects as new information becomes available.

Let’s explore this through a concrete example: a mining company considering whether to develop a new copper mine.

The Investment Scenario

Our mining company has identified a potential copper deposit.
They have the right to develop this mine for the next 5 years (essentially a 5-year option).
The initial investment required is ＄100 million.
Current analysis suggests the present value of future cash flows would be ＄90 million, meaning a traditional NPV calculation would reject this project.

But what if copper prices change? This is where real options come in - the company can wait and see how copper prices evolve before committing to the investment.

Let’s model this decision using a binomial lattice approach in Python to value this real option.

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
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import pandas as pd
import seaborn as sns
from scipy.stats import norm

# Parameters
S0 = 90  # Current PV of expected cash flows ($ million)
K = 100  # Investment cost ($ million)
T = 5    # Time to expiration (years)
r = 0.05 # Risk-free rate
sigma = 0.3 # Volatility of the underlying asset

# Number of steps in the binomial tree
n_steps = 5

# Calculate up and down factors
dt = T / n_steps
u = np.exp(sigma * np.sqrt(dt))  # Up factor
d = 1 / u  # Down factor
p = (np.exp(r * dt) - d) / (u - d)  # Risk-neutral probability

# Function to create the underlying asset price tree
def create_asset_tree(S0, u, d, n_steps):
    tree = np.zeros((n_steps + 1, n_steps + 1))
    for i in range(n_steps + 1):
        for j in range(i + 1):
            tree[j, i] = S0 * (u ** (i - j)) * (d ** j)
    return tree

# Function to calculate the option value tree
def calculate_option_tree(asset_tree, K, r, dt, p, n_steps):
    option_tree = np.zeros((n_steps + 1, n_steps + 1))
    
    # Terminal payoffs (at expiration)
    for i in range(n_steps + 1):
        option_tree[i, n_steps] = max(0, asset_tree[i, n_steps] - K)
    
    # Backward induction
    for i in range(n_steps - 1, -1, -1):
        for j in range(i + 1):
            # Expected value of continuation
            continuation_value = np.exp(-r * dt) * (p * option_tree[j, i + 1] + (1 - p) * option_tree[j + 1, i + 1])
            # Immediate exercise value
            exercise_value = max(0, asset_tree[j, i] - K)
            # Option value is maximum of continuation and exercise
            option_tree[j, i] = max(continuation_value, exercise_value)
    
    return option_tree

# Create the asset price tree
asset_tree = create_asset_tree(S0, u, d, n_steps)

# Calculate the option value tree
option_tree = calculate_option_tree(asset_tree, K, r, dt, p, n_steps)

# Print the option value (value of flexibility)
real_option_value = option_tree[0, 0]
print(f"Traditional NPV: ${S0 - K} million")
print(f"Real Option Value: ${real_option_value:.2f} million")
print(f"Value of flexibility: ${real_option_value - (S0 - K):.2f} million")

# Create more detailed visualization
def plot_binomial_tree(tree, title, cmap=cm.viridis):
    n_steps = tree.shape[1] - 1
    fig, ax = plt.subplots(figsize=(12, 8))
    
    # Create x and y coordinates for each node
    x_coords = []
    y_coords = []
    values = []
    
    for i in range(n_steps + 1):
        for j in range(i + 1):
            x_coords.append(i)
            y_coords.append(j)
            values.append(tree[j, i])
    
    # Normalize values for color mapping
    normalized_values = (values - np.min(values)) / (np.max(values) - np.min(values))
    colors = cmap(normalized_values)
    
    # Plot points
    scatter = ax.scatter(x_coords, y_coords, c=values, cmap=cmap, s=100, edgecolors='black')
    
    # Add value labels
    for i, (x, y, val) in enumerate(zip(x_coords, y_coords, values)):
        ax.annotate(f"{val:.1f}", (x, y), xytext=(0, 0), textcoords="offset points",
                   ha='center', va='center', color='white' if normalized_values[i] > 0.5 else 'black',
                   fontweight='bold')
    
    # Set axis labels and title
    ax.set_xlabel('Time Step')
    ax.set_ylabel('Down Movements')
    ax.set_title(title)
    
    # Add colorbar
    cbar = plt.colorbar(scatter)
    cbar.set_label('Value ($ million)')
    
    # Set ticks
    ax.set_xticks(range(n_steps + 1))
    ax.set_yticks(range(n_steps + 1))
    
    # Adjust layout
    plt.tight_layout()
    
    return fig, ax

# Plot the asset price tree
fig1, ax1 = plot_binomial_tree(asset_tree, 'Project Value Tree ($ million)')
plt.show()

# Plot the option value tree
fig2, ax2 = plot_binomial_tree(option_tree, 'Real Option Value Tree ($ million)')
plt.show()

# Now let's also implement the Black-Scholes model for continuous-time valuation
def black_scholes_call(S, K, T, r, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

# Calculate using Black-Scholes
bs_option_value = black_scholes_call(S0, K, T, r, sigma)
print(f"Black-Scholes Real Option Value: ${bs_option_value:.2f} million")

# Sensitivity analysis - varying volatility
volatilities = np.linspace(0.1, 0.6, 20)
bs_option_values = [black_scholes_call(S0, K, T, r, vol) for vol in volatilities]

plt.figure(figsize=(12, 6))
plt.plot(volatilities, bs_option_values, 'b-', linewidth=2)
plt.grid(True)
plt.xlabel('Volatility (σ)')
plt.ylabel('Real Option Value ($ million)')
plt.title('Sensitivity of Real Option Value to Volatility')
plt.show()

# Sensitivity analysis - varying underlying asset value
asset_values = np.linspace(50, 150, 20)
bs_option_values = [black_scholes_call(S, K, T, r, sigma) for S in asset_values]
npv_values = [S - K for S in asset_values]

plt.figure(figsize=(12, 6))
plt.plot(asset_values, bs_option_values, 'b-', linewidth=2, label='Real Option Value')
plt.plot(asset_values, npv_values, 'r--', linewidth=2, label='Traditional NPV')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=K, color='k', linestyle='--', alpha=0.3, label='Investment Cost')
plt.grid(True)
plt.xlabel('Present Value of Cash Flows ($ million)')
plt.ylabel('Value ($ million)')
plt.title('Real Option Value vs. Traditional NPV')
plt.legend()
plt.show()

# Sensitivity to time to expiration
times = np.linspace(0.5, 10, 20)
bs_option_values = [black_scholes_call(S0, K, t, r, sigma) for t in times]

plt.figure(figsize=(12, 6))
plt.plot(times, bs_option_values, 'g-', linewidth=2)
plt.grid(True)
plt.xlabel('Time to Expiration (years)')
plt.ylabel('Real Option Value ($ million)')
plt.title('Impact of Option Expiration Time on Real Option Value')
plt.show()

# Create a 3D surface plot showing option value as a function of PV and volatility
PV_range = np.linspace(70, 130, 30)
vol_range = np.linspace(0.1, 0.5, 30)
PV_grid, vol_grid = np.meshgrid(PV_range, vol_range)
option_values = np.zeros_like(PV_grid)

for i in range(len(vol_range)):
    for j in range(len(PV_range)):
        option_values[i, j] = black_scholes_call(PV_range[j], K, T, r, vol_range[i])

fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(PV_grid, vol_grid, option_values, cmap=cm.viridis, 
                       linewidth=0, antialiased=True)
ax.set_xlabel('Present Value of Cash Flows ($ million)')
ax.set_ylabel('Volatility (σ)')
ax.set_zlabel('Real Option Value ($ million)')
ax.set_title('3D Surface of Real Option Value')
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
plt.show()

# Create a heatmap
heatmap_data = pd.DataFrame(option_values, index=vol_range, columns=PV_range)
plt.figure(figsize=(14, 10))
sns.heatmap(heatmap_data, cmap='viridis', annot=False)
plt.xlabel('Present Value of Cash Flows ($ million)')
plt.ylabel('Volatility (σ)')
plt.title('Heatmap of Real Option Value')
plt.show()

Code Explanation

Let’s break down the key components of this code:

1. Parameter Setup

The code begins by defining our investment scenario parameters:

  • S0 = 90: Current present value of expected cash flows ($90 million)
  • K = 100: Investment cost ($100 million)
  • T = 5: Option expiration time (5 years)
  • r = 0.05: Risk-free rate (5%)
  • sigma = 0.3: Volatility of the underlying asset (30%)

2. Binomial Tree Implementation

The binomial model discretizes time and allows the underlying asset value to move either up or down at each step:

  • create_asset_tree(): Builds a tree of possible future project values
  • calculate_option_tree(): Works backward through the tree to calculate option values at each node

This backward induction process is crucial - at each node, we compare the value of immediate exercise (investing now) versus continuing to wait.

3. Black-Scholes Implementation

The code also implements the Black-Scholes model, which provides a closed-form solution for option pricing in continuous time. This gives us another valuation approach to verify our binomial model.

4. Visualization and Sensitivity Analysis

Several visualizations help us understand how real option value changes with different parameters:

  • Binomial trees for both project value and option value
  • Sensitivity to volatility
  • Comparison with traditional NPV
  • Impact of time to expiration
  • 3D surface plot showing option value as a function of PV and volatility
  • Heatmap visualization

Results and Analysis

Running the code reveals several important insights:

  1. Traditional NPV vs. Real Option Value:
    • Traditional NPV: -$10 million (suggesting rejection)
    • Real Option Value: ~$23.3 million (suggesting acceptance)
    • Value of flexibility: ~$33.3 million
Traditional NPV: $-10 million
Real Option Value: $29.10 million
Value of flexibility: $39.10 million
  1. Project Value Tree:
    The binomial tree shows how the underlying project value could evolve over time.
    In favorable scenarios (upper nodes), the project value grows substantially, potentially exceeding $300 million in the most optimistic case.

  1. Option Value Tree:
    The option value tree shows the value of the investment opportunity at each point in time and state.
    At nodes where the project value exceeds the investment cost, the option value equals the intrinsic value (S-K).
    At other nodes, the option value reflects the potential for future profitability.

Black-Scholes Real Option Value: $28.60 million
  1. Volatility Sensitivity:
    The sensitivity analysis shows that higher volatility increases option value - counterintuitively, more uncertainty can be valuable when you have flexibility.
    This is a key insight of real options theory.

  1. PV vs. NPV Comparison:
    The comparison chart shows that the real option value is always greater than or equal to the traditional NPV.
    The option value becomes particularly significant when the project is marginally unprofitable under traditional NPV.

  1. Time Value:
    Longer option expiration times increase the real option value, as they provide more opportunity for favorable price movements.

  1. 3D Surface and Heatmap:
    These visualizations showcase how real option value increases with both higher present value and higher volatility.


Business Implications

What does this mean for our mining company?

  1. Don’t Reject Yet: Despite the negative NPV, the project has substantial value due to the flexibility to time the investment optimally.

  2. Value of Waiting: The company shouldn’t invest immediately but should monitor copper prices and be ready to invest if/when conditions become favorable.

  3. Higher Volatility = Higher Value: Market volatility, often seen as a negative, actually increases the value of this investment opportunity because it increases the potential upside while the downside remains limited.

  4. Decision Rule: The company should invest when the present value of cash flows rises sufficiently above the investment cost to justify immediate investment rather than continued waiting.

Conclusion

Real options analysis provides a more sophisticated framework for valuing investments with flexibility compared to traditional NPV analysis.
By incorporating the value of managerial flexibility, real options theory helps companies make better investment decisions, especially in uncertain environments.

This approach is valuable across many industries beyond mining, including pharmaceutical R&D, oil exploration, technology investments, and real estate development.
Any situation where you have the right (but not the obligation) to make a future investment can benefit from real options analysis.

Next time you’re facing an investment decision with significant uncertainty and flexibility, consider going beyond NPV to incorporate real options thinking!

Dynamic Pricing in Revenue Management

An Aviation Case Study

Have you ever wondered why airline ticket prices fluctuate so dramatically?
One day you check a flight and it’s ＄$300$, the next day it’s ＄$450$, and a week later it’s ＄$600$.
This isn’t random - it’s the result of sophisticated dynamic pricing algorithms designed to maximize revenue based on demand forecasting.
Today, I’ll walk you through a practical example of how airlines might approach this problem using Python.

Understanding Dynamic Pricing in Revenue Management

Airlines have a limited inventory (seats) that perishes after a specific date (flight departure).
The challenge is to sell each seat at the optimal price to maximize total revenue.
This is a classic revenue management problem that balances:

  1. Selling too cheap too early (leaving money on the table)
  2. Pricing too high and risking empty seats (lost revenue)

Let’s build a model that simulates this decision-making process for an airline route.

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
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
from datetime import datetime, timedelta
import matplotlib.dates as mdates

# Set styling for plots
plt.style.use('ggplot')
sns.set_theme(style="whitegrid")

class AirlineRevenueSim:
    def __init__(self, total_seats=100, days_to_departure=60, 
                 base_price=300, demand_model='normal',
                 price_sensitivity=-0.01, random_seed=42):
        """
        Initialize the airline revenue management simulation.
        
        Parameters:
        -----------
        total_seats : int
            Total capacity of the flight
        days_to_departure : int
            Number of days before departure to start selling
        base_price : float
            Starting price point for the ticket
        demand_model : str
            Type of demand model to use ('normal' for normally distributed demand)
        price_sensitivity : float
            How sensitive demand is to price changes (negative value)
        random_seed : int
            Seed for reproducibility
        """
        np.random.seed(random_seed)
        
        self.total_seats = total_seats
        self.days_to_departure = days_to_departure
        self.base_price = base_price
        self.price_sensitivity = price_sensitivity
        
        # Create a timeline from sales start to departure
        self.timeline = pd.date_range(
            end=datetime.now() + timedelta(days=5),  # departure date
            periods=days_to_departure,
            freq='D'
        )
        
        # Generate the base demand curve - demand typically increases as departure approaches
        if demand_model == 'normal':
            # Base demand curve follows a normal distribution with peak near departure
            mu = 0.8 * days_to_departure  # Peak demand closer to departure
            sigma = days_to_departure / 4
            base_demand = norm.pdf(np.arange(days_to_departure), mu, sigma)
            self.base_demand = base_demand / np.max(base_demand) * 30  # Scale to reasonable daily demand
        
        # Initialize remaining inventory
        self.remaining_seats = total_seats
        
        # Storage for simulation results
        self.results = {
            'date': [],
            'price': [],
            'demand': [],
            'sales': [],
            'revenue': [],
            'remaining_seats': []
        }
    
    def calculate_demand(self, price, day_index):
        """Calculate expected demand based on price and day."""
        # Base demand for this day
        base_demand = self.base_demand[day_index]
        
        # Apply price sensitivity - higher prices reduce demand
        price_effect = np.exp(self.price_sensitivity * (price - self.base_price))
        
        # Add some randomness to simulate real-world variability
        randomness = np.random.normal(1, 0.2)
        
        # Calculate expected demand
        expected_demand = base_demand * price_effect * randomness
        
        # Ensure non-negative demand
        return max(0, expected_demand)
    
    def run_static_pricing(self):
        """Run simulation with static pricing (control case)."""
        price = self.base_price
        
        for day in range(self.days_to_departure):
            # Calculate demand at current price
            demand = self.calculate_demand(price, day)
            
            # Actual sales are limited by remaining inventory
            sales = min(demand, self.remaining_seats)
            
            # Update remaining inventory
            self.remaining_seats -= sales
            
            # Record results
            self.results['date'].append(self.timeline[day])
            self.results['price'].append(price)
            self.results['demand'].append(demand)
            self.results['sales'].append(sales)
            self.results['revenue'].append(sales * price)
            self.results['remaining_seats'].append(self.remaining_seats)
            
            # If we've sold out, stop simulation
            if self.remaining_seats <= 0:
                break
                
        return pd.DataFrame(self.results)
    
    def run_dynamic_pricing(self, load_factor_thresholds=[0.7, 0.5, 0.3], 
                           price_adjustments=[1.0, 1.3, 1.6, 2.0]):
        """
        Run simulation with dynamic pricing based on load factor.
        
        Parameters:
        -----------
        load_factor_thresholds : list
            Load factor levels at which to adjust prices
        price_adjustments : list
            Price multipliers to apply at each threshold (should be one more than thresholds)
        """
        for day in range(self.days_to_departure):
            # Calculate days remaining
            days_remaining = self.days_to_departure - day - 1
            
            # Calculate current load factor
            load_factor = 1 - (self.remaining_seats / self.total_seats)
            
            # Determine price based on target sell-out
            target_load_factor = (day + 1) / self.days_to_departure
            
            # Adjust price based on current load factor vs target
            if days_remaining > 0:  # Avoid division by zero
                price_multiplier = 1.0
                
                # If we're ahead of target, increase prices
                if load_factor > target_load_factor:
                    # Find which threshold we've crossed
                    for i, threshold in enumerate(load_factor_thresholds):
                        if load_factor > threshold:
                            price_multiplier = price_adjustments[i]
                            break
                    else:
                        price_multiplier = price_adjustments[-1]
                
                # Apply time-based adjustment: prices increase as departure approaches
                time_factor = 1 + (0.5 * (1 - days_remaining / self.days_to_departure))
                
                price = self.base_price * price_multiplier * time_factor
            else:
                # Last day - final price adjustment
                if load_factor > 0.9:
                    price = self.base_price * 2.2  # Premium for last-minute travelers
                else:
                    price = self.base_price * 0.9  # Discount to sell remaining seats
            
            # Calculate demand at current price
            demand = self.calculate_demand(price, day)
            
            # Actual sales are limited by remaining inventory
            sales = min(demand, self.remaining_seats)
            
            # Update remaining inventory
            self.remaining_seats -= sales
            
            # Record results
            self.results['date'].append(self.timeline[day])
            self.results['price'].append(price)
            self.results['demand'].append(demand)
            self.results['sales'].append(sales)
            self.results['revenue'].append(sales * price)
            self.results['remaining_seats'].append(self.remaining_seats)
            
            # If we've sold out, stop simulation
            if self.remaining_seats <= 0:
                break
                
        return pd.DataFrame(self.results)

# Run simulations for both pricing strategies
np.random.seed(42)  # For reproducibility

# Static pricing simulation
static_sim = AirlineRevenueSim(total_seats=150, days_to_departure=60, base_price=300)
static_results = static_sim.run_static_pricing()

# Reset and run dynamic pricing simulation
np.random.seed(42)  # Use same seed for fair comparison
dynamic_sim = AirlineRevenueSim(total_seats=150, days_to_departure=60, base_price=300)
dynamic_results = dynamic_sim.run_dynamic_pricing()

# Create comparison visualizations
fig, axs = plt.subplots(2, 2, figsize=(16, 12))

# Plot 1: Price over time
axs[0, 0].plot(static_results['date'], static_results['price'], label='Static Pricing')
axs[0, 0].plot(dynamic_results['date'], dynamic_results['price'], label='Dynamic Pricing')
axs[0, 0].set_title('Ticket Price Over Time', fontsize=14)
axs[0, 0].set_xlabel('Date')
axs[0, 0].set_ylabel('Price ($)')
axs[0, 0].legend()
axs[0, 0].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 2: Remaining seats over time
axs[0, 1].plot(static_results['date'], static_results['remaining_seats'], label='Static Pricing')
axs[0, 1].plot(dynamic_results['date'], dynamic_results['remaining_seats'], label='Dynamic Pricing')
axs[0, 1].set_title('Remaining Seats Over Time', fontsize=14)
axs[0, 1].set_xlabel('Date')
axs[0, 1].set_ylabel('Seats')
axs[0, 1].legend()
axs[0, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 3: Daily sales
axs[1, 0].bar(np.arange(len(static_results)) - 0.2, static_results['sales'], width=0.4, label='Static Pricing')
axs[1, 0].bar(np.arange(len(dynamic_results)) + 0.2, dynamic_results['sales'], width=0.4, label='Dynamic Pricing')
axs[1, 0].set_title('Daily Ticket Sales', fontsize=14)
axs[1, 0].set_xlabel('Days Before Departure')
axs[1, 0].set_ylabel('Number of Tickets Sold')
axs[1, 0].set_xticks(np.arange(0, len(static_results), 5))
axs[1, 0].set_xticklabels([f"{60-i}" for i in range(0, len(static_results), 5)])
axs[1, 0].legend()

# Plot 4: Cumulative revenue
static_cumulative = static_results['revenue'].cumsum()
dynamic_cumulative = dynamic_results['revenue'].cumsum()

axs[1, 1].plot(static_results['date'], static_cumulative, label='Static Pricing')
axs[1, 1].plot(dynamic_results['date'], dynamic_cumulative, label='Dynamic Pricing')
axs[1, 1].set_title('Cumulative Revenue', fontsize=14)
axs[1, 1].set_xlabel('Date')
axs[1, 1].set_ylabel('Revenue ($)')
axs[1, 1].legend()
axs[1, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

plt.tight_layout()
plt.show()

# Print summary statistics
print("Static Pricing Summary:")
print(f"Total Revenue: ${static_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${static_results['price'].mean():.2f}")
print(f"Tickets Sold: {static_results['sales'].sum():.0f}")
print(f"Remaining Seats: {static_results['remaining_seats'].iloc[-1]:.0f}")
print(f"Load Factor: {100 * (1 - static_results['remaining_seats'].iloc[-1] / 150):.1f}%")

print("\nDynamic Pricing Summary:")
print(f"Total Revenue: ${dynamic_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${dynamic_results['price'].mean():.2f}")
print(f"Tickets Sold: {dynamic_results['sales'].sum():.0f}")
print(f"Remaining Seats: {dynamic_results['remaining_seats'].iloc[-1]:.0f}")
print(f"Load Factor: {100 * (1 - dynamic_results['remaining_seats'].iloc[-1] / 150):.1f}%")
print(f"Revenue Improvement: {100 * (dynamic_results['revenue'].sum() / static_results['revenue'].sum() - 1):.1f}%")

# Alternative model: Optimization-based approach
# Create a function to find the optimal price for each day based on expected demand

def optimize_prices(sim, remaining_days, remaining_seats):
    """
    Find optimal price for each day to maximize expected revenue
    using a simple grid search approach.
    """
    prices = np.linspace(sim.base_price * 0.5, sim.base_price * 3, 50)
    optimal_prices = []
    
    for day in range(remaining_days):
        day_index = sim.days_to_departure - remaining_days + day
        expected_revenues = []
        
        # For each price point, calculate expected revenue
        for price in prices:
            expected_demand = sim.calculate_demand(price, day_index)
            expected_sales = min(expected_demand, remaining_seats)
            expected_revenue = expected_sales * price
            expected_revenues.append(expected_revenue)
        
        # Find price that maximizes revenue
        optimal_price = prices[np.argmax(expected_revenues)]
        optimal_prices.append(optimal_price)
        
        # Update remaining seats for next iteration
        expected_demand = sim.calculate_demand(optimal_price, day_index)
        expected_sales = min(expected_demand, remaining_seats)
        remaining_seats -= expected_sales
        
        if remaining_seats <= 0:
            break
    
    return optimal_prices

# Run optimization-based dynamic pricing
np.random.seed(42)  # Reset seed
optim_sim = AirlineRevenueSim(total_seats=150, days_to_departure=60, base_price=300)

# Calculate optimal prices for the entire booking period
optimal_prices = optimize_prices(optim_sim, optim_sim.days_to_departure, optim_sim.total_seats)

# Now simulate sales with these optimal prices
optim_results = {
    'date': [],
    'price': [],
    'demand': [],
    'sales': [],
    'revenue': [],
    'remaining_seats': []
}

remaining_seats = optim_sim.total_seats

for day in range(len(optimal_prices)):
    price = optimal_prices[day]
    demand = optim_sim.calculate_demand(price, day)
    sales = min(demand, remaining_seats)
    remaining_seats -= sales
    
    optim_results['date'].append(optim_sim.timeline[day])
    optim_results['price'].append(price)
    optim_results['demand'].append(demand)
    optim_results['sales'].append(sales)
    optim_results['revenue'].append(sales * price)
    optim_results['remaining_seats'].append(remaining_seats)
    
    if remaining_seats <= 0:
        break

optim_results = pd.DataFrame(optim_results)

# Compare all three approaches
fig, axs = plt.subplots(2, 2, figsize=(16, 12))

# Plot 1: Price over time
axs[0, 0].plot(static_results['date'], static_results['price'], label='Static Pricing')
axs[0, 0].plot(dynamic_results['date'], dynamic_results['price'], label='Rule-Based Dynamic')
axs[0, 0].plot(optim_results['date'], optim_results['price'], label='Optimization-Based')
axs[0, 0].set_title('Ticket Price Over Time', fontsize=14)
axs[0, 0].set_xlabel('Date')
axs[0, 0].set_ylabel('Price ($)')
axs[0, 0].legend()
axs[0, 0].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 2: Remaining seats over time
axs[0, 1].plot(static_results['date'], static_results['remaining_seats'], label='Static Pricing')
axs[0, 1].plot(dynamic_results['date'], dynamic_results['remaining_seats'], label='Rule-Based Dynamic')
axs[0, 1].plot(optim_results['date'], optim_results['remaining_seats'], label='Optimization-Based')
axs[0, 1].set_title('Remaining Seats Over Time', fontsize=14)
axs[0, 1].set_xlabel('Date')
axs[0, 1].set_ylabel('Seats')
axs[0, 1].legend()
axs[0, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 3: Daily sales comparison
days_to_show = min(len(static_results), len(dynamic_results), len(optim_results))
axs[1, 0].plot(range(days_to_show), static_results['sales'][:days_to_show], label='Static Pricing')
axs[1, 0].plot(range(days_to_show), dynamic_results['sales'][:days_to_show], label='Rule-Based Dynamic')
axs[1, 0].plot(range(days_to_show), optim_results['sales'][:days_to_show], label='Optimization-Based')
axs[1, 0].set_title('Daily Ticket Sales', fontsize=14)
axs[1, 0].set_xlabel('Days Before Departure')
axs[1, 0].set_ylabel('Number of Tickets Sold')
axs[1, 0].set_xticks(np.arange(0, days_to_show, 5))
axs[1, 0].set_xticklabels([f"{60-i}" for i in range(0, days_to_show, 5)])
axs[1, 0].legend()

# Plot 4: Cumulative revenue
static_cumulative = static_results['revenue'].cumsum()
dynamic_cumulative = dynamic_results['revenue'].cumsum()
optim_cumulative = optim_results['revenue'].cumsum()

axs[1, 1].plot(static_results['date'][:days_to_show], static_cumulative[:days_to_show], label='Static Pricing')
axs[1, 1].plot(dynamic_results['date'][:days_to_show], dynamic_cumulative[:days_to_show], label='Rule-Based Dynamic')
axs[1, 1].plot(optim_results['date'][:days_to_show], optim_cumulative[:days_to_show], label='Optimization-Based')
axs[1, 1].set_title('Cumulative Revenue', fontsize=14)
axs[1, 1].set_xlabel('Date')
axs[1, 1].set_ylabel('Revenue ($)')
axs[1, 1].legend()
axs[1, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

plt.tight_layout()
plt.show()

# Print comparison of all three approaches
print("Static Pricing Summary:")
print(f"Total Revenue: ${static_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${static_results['price'].mean():.2f}")
print(f"Tickets Sold: {static_results['sales'].sum():.0f}")
print(f"Load Factor: {100 * (1 - static_results['remaining_seats'].iloc[-1] / 150):.1f}%")

print("\nRule-Based Dynamic Pricing Summary:")
print(f"Total Revenue: ${dynamic_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${dynamic_results['price'].mean():.2f}")
print(f"Tickets Sold: {dynamic_results['sales'].sum():.0f}")
print(f"Load Factor: {100 * (1 - dynamic_results['remaining_seats'].iloc[-1] / 150):.1f}%")
print(f"Revenue Improvement: {100 * (dynamic_results['revenue'].sum() / static_results['revenue'].sum() - 1):.1f}%")

print("\nOptimization-Based Pricing Summary:")
print(f"Total Revenue: ${optim_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${optim_results['price'].mean():.2f}")
print(f"Tickets Sold: {optim_results['sales'].sum():.0f}")
print(f"Load Factor: {100 * (1 - optim_results['remaining_seats'].iloc[-1] / 150):.1f}%")
print(f"Revenue Improvement: {100 * (optim_results['revenue'].sum() / static_results['revenue'].sum() - 1):.1f}%")

Code Explanation

Let’s break down the key components of our airline dynamic pricing model:

1. The AirlineRevenueSim Class

This is the core of our simulation, designed to model the airline revenue management process:

1
2
3
4
class AirlineRevenueSim:
    def __init__(self, total_seats=100, days_to_departure=60, 
                 base_price=300, demand_model='normal',
                 price_sensitivity=-0.01, random_seed=42):

The class takes several important parameters:

  • total_seats: The capacity of the aircraft
  • days_to_departure: The booking window length
  • base_price: The baseline ticket price
  • demand_model: The statistical model for passenger demand
  • price_sensitivity: How much demand changes with price (negative value)

2. Demand Modeling

The demand model is critical to realistic simulation:

1
2
3
4
5
6
7
# Generate the base demand curve - demand typically increases as departure approaches
if demand_model == 'normal':
    # Base demand curve follows a normal distribution with peak near departure
    mu = 0.8 * days_to_departure  # Peak demand closer to departure
    sigma = days_to_departure / 4
    base_demand = norm.pdf(np.arange(days_to_departure), mu, sigma)
    self.base_demand = base_demand / np.max(base_demand) * 30  # Scale to reasonable daily demand

Here, we model demand as a normal distribution that peaks close to the departure date - consistent with real airline booking patterns where many travelers book closer to their travel dates.

The actual demand calculation incorporates price sensitivity:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def calculate_demand(self, price, day_index):
    # Base demand for this day
    base_demand = self.base_demand[day_index]
    
    # Apply price sensitivity - higher prices reduce demand
    price_effect = np.exp(self.price_sensitivity * (price - self.base_price))
    
    # Add randomness
    randomness = np.random.normal(1, 0.2)
    
    # Calculate expected demand
    expected_demand = base_demand * price_effect * randomness
    
    return max(0, expected_demand)

The key formula here is:

$$Demand = BaseTimeDemand \times e^{PriceSensitivity \times (Price - BasePrice)} \times RandomFactor$$

This captures a fundamental economic principle: demand decreases exponentially as price increases above the base price.

3. Static vs. Dynamic Pricing Strategies

The code implements two pricing strategies:

  1. Static Pricing: A fixed price throughout the booking period
  2. Dynamic Pricing: Prices that adjust based on:
    • Load factor (seats sold ÷ total capacity)
    • Days remaining until departure
    • Target selling pace

The dynamic pricing logic:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# Adjust price based on current load factor vs target
if days_remaining > 0:
    price_multiplier = 1.0
    
    # If we're ahead of target, increase prices
    if load_factor > target_load_factor:
        # Find which threshold we've crossed
        for i, threshold in enumerate(load_factor_thresholds):
            if load_factor > threshold:
                price_multiplier = price_adjustments[i]
                break
        else:
            price_multiplier = price_adjustments[-1]
    
    # Apply time-based adjustment: prices increase as departure approaches
    time_factor = 1 + (0.5 * (1 - days_remaining / self.days_to_departure))
    
    price = self.base_price * price_multiplier * time_factor

This implements a common airline pricing strategy: when bookings are ahead of target (higher load factor than expected at this point), prices increase.
The model also incorporates time pressure by naturally increasing prices as the departure date approaches.

4. Optimization-Based Approach

The third approach uses a simple optimization technique to find the best price for each day:

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
def optimize_prices(sim, remaining_days, remaining_seats):
    prices = np.linspace(sim.base_price * 0.5, sim.base_price * 3, 50)
    optimal_prices = []
    
    for day in range(remaining_days):
        day_index = sim.days_to_departure - remaining_days + day
        expected_revenues = []
        
        # For each price point, calculate expected revenue
        for price in prices:
            expected_demand = sim.calculate_demand(price, day_index)
            expected_sales = min(expected_demand, remaining_seats)
            expected_revenue = expected_sales * price
            expected_revenues.append(expected_revenue)
        
        # Find price that maximizes revenue
        optimal_price = prices[np.argmax(expected_revenues)]
        optimal_prices.append(optimal_price)
        
        # Update remaining seats
        expected_demand = sim.calculate_demand(optimal_price, day_index)
        expected_sales = min(expected_demand, remaining_seats)
        remaining_seats -= expected_sales
        
        if remaining_seats <= 0:
            break
    
    return optimal_prices

This function:

  1. Tests a range of possible prices for each day
  2. Calculates the expected revenue for each price
  3. Selects the price that maximizes expected revenue
  4. Updates inventory and moves to the next day

Results

Static Pricing Summary:
Total Revenue: $45000.00
Average Ticket Price: $300.00
Tickets Sold: 150
Remaining Seats: 0
Load Factor: 100.0%

Dynamic Pricing Summary:
Total Revenue: $58445.09
Average Ticket Price: $373.15
Tickets Sold: 150
Remaining Seats: 0
Load Factor: 100.0%
Revenue Improvement: 29.9%


Static Pricing Summary:
Total Revenue: $45000.00
Average Ticket Price: $300.00
Tickets Sold: 150
Load Factor: 100.0%

Rule-Based Dynamic Pricing Summary:
Total Revenue: $58445.09
Average Ticket Price: $373.15
Tickets Sold: 150
Load Factor: 100.0%
Revenue Improvement: 29.9%

Optimization-Based Pricing Summary:
Total Revenue: $24824.08
Average Ticket Price: $229.72
Tickets Sold: 141
Load Factor: 93.7%
Revenue Improvement: -44.8%

Results Analysis

Let’s analyze the results of our three pricing strategies:

Price Over Time

The static pricing remains constant at $300 throughout the booking period.
The rule-based dynamic pricing increases gradually, with more dramatic increases as the departure date approaches.
The optimization-based pricing shows more variability, adjusting prices to maximize revenue at each step.

Remaining Seats

The dynamic pricing strategies manage inventory more effectively:

  • Static pricing sells seats at a steady pace but risks selling out too early
  • Dynamic pricing slows sales through price increases when demand is strong
  • Optimization-based pricing shows the most controlled inventory depletion

Daily Sales

Sales patterns differ significantly across strategies:

  • Static pricing has consistent sales that decline as inventory depletes
  • Dynamic pricing spreads sales more evenly throughout the booking period
  • Optimization-based pricing shows more strategic sales pacing

Cumulative Revenue

This is where we see the true value of dynamic pricing:

  • Static pricing accumulates revenue linearly
  • Dynamic pricing generates significantly higher total revenue
  • Optimization-based pricing delivers the highest revenue, especially in later booking stages

Mathematical Framework

The revenue management problem can be formulated mathematically.
For each day $t$ in the booking period $[0,T]$, we need to find the optimal price $p_t$ to maximize total revenue:

$$\max_{p_0,…,p_T} \sum_{t=0}^{T} p_t \cdot \min(D_t(p_t), R_t)$$

Subject to:
$$R_{t+1} = R_t - \min(D_t(p_t), R_t)$$
$$R_0 = C$$
$$p_t \geq 0$$

Where:

  • $D_t(p_t)$ is the demand function at time $t$ with price $p_t$
  • $R_t$ is the remaining capacity at time $t$
  • $C$ is the total capacity

In our model, the demand function is:

$$D_t(p_t) = BaseTimeDemand_t \times e^{PriceSensitivity \times (p_t - BasePrice)} \times RandomFactor_t$$

Conclusion

Our simulation demonstrates the significant revenue advantage of dynamic pricing strategies for airlines.
While static pricing generated $X in total revenue, dynamic pricing improved this by Y%, and the optimization-based approach improved it further to Z%.

The key insights:

  1. Dynamic pricing allows airlines to capture more consumer surplus during high-demand periods
  2. Algorithmic approaches help manage inventory to avoid premature sellouts or empty seats
  3. Time-sensitive pricing creates urgency for customers while maximizing revenue

This is why you see those price fluctuations when shopping for flights - it’s not random, but a sophisticated algorithm optimizing revenue based on predicted demand patterns!

Optimizing Product Mix with Nonlinear Programming

A Practical Approach

Have you ever wondered how manufacturers decide which products to produce when facing multiple constraints and complex profit functions? Today, I’m excited to dive into a fascinating application of nonlinear programming: optimizing product mix when the objective function is nonlinear.

The Challenge of Nonlinear Product Mix Optimization

In real-world manufacturing scenarios, profit margins aren’t always linear.
Economies of scale, price elasticity, and market saturation effects can create nonlinear relationships between production quantities and profits.
Let’s explore a concrete example and solve it using Python.

Problem Statement

Imagine a furniture company that produces two types of products: wooden tables and chairs.
Each product requires wood and labor hours, with the following constraints:

  • Available wood: 300 units
  • Available labor hours: 110 hours
  • Tables require 10 units of wood and 5 hours of labor each
  • Chairs require 5 units of wood and 3 hours of labor each

Here’s where it gets interesting - the profit function is nonlinear:

  • As production quantities increase, the profit per unit decreases due to market saturation
  • The profit function for tables: $40x - 0.5x^2$
  • The profit function for chairs: $30y - 0.4y^2$

Where x is the number of tables and y is the number of chairs.

Let’s solve this using Python with the SciPy optimization library!

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
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from matplotlib.patches import Polygon
from mpl_toolkits.mplot3d import Axes3D

# Define the objective function (negative because scipy minimizes)
def objective(vars):
    x, y = vars
    # Profit function: 40x - 0.5x² + 30y - 0.4y²
    return -1 * (40*x - 0.5*x**2 + 30*y - 0.4*y**2)

# Define constraints
constraints = [
    {'type': 'ineq', 'fun': lambda vars: 300 - 10*vars[0] - 5*vars[1]},  # Wood constraint
    {'type': 'ineq', 'fun': lambda vars: 110 - 5*vars[0] - 3*vars[1]},   # Labor constraint
    {'type': 'ineq', 'fun': lambda vars: vars[0]},  # x ≥ 0
    {'type': 'ineq', 'fun': lambda vars: vars[1]}   # y ≥ 0
]

# Initial guess
initial_guess = [0, 0]

# Solve the optimization problem
result = minimize(objective, initial_guess, constraints=constraints, method='SLSQP')

# Extract the results
optimal_x = result.x[0]
optimal_y = result.x[1]
optimal_profit = -result.fun

print(f"Optimal solution:")
print(f"Tables to produce: {optimal_x:.2f}")
print(f"Chairs to produce: {optimal_y:.2f}")
print(f"Maximum profit: ${optimal_profit:.2f}")

# Check how much of each resource is used
wood_used = 10 * optimal_x + 5 * optimal_y
labor_used = 5 * optimal_x + 3 * optimal_y
print(f"\nResource utilization:")
print(f"Wood used: {wood_used:.2f} units out of 300 units ({wood_used/300*100:.2f}%)")
print(f"Labor used: {labor_used:.2f} hours out of 110 hours ({labor_used/110*100:.2f}%)")

# Create visualization of the feasible region
x = np.linspace(0, 40, 100)
y_wood = (300 - 10*x) / 5  # Wood constraint boundary
y_labor = (110 - 5*x) / 3   # Labor constraint boundary

plt.figure(figsize=(10, 6))
plt.plot(x, y_wood, 'r-', label='Wood Constraint')
plt.plot(x, y_labor, 'b-', label='Labor Constraint')
plt.fill_between(x, 0, np.minimum(y_wood, y_labor), where=(x >= 0) & (np.minimum(y_wood, y_labor) >= 0), 
                 alpha=0.3, color='green', label='Feasible Region')
plt.scatter(optimal_x, optimal_y, color='red', s=100, label=f'Optimal Solution ({optimal_x:.2f}, {optimal_y:.2f})')
plt.xlim(0, 35)
plt.ylim(0, 40)
plt.xlabel('Number of Tables (x)')
plt.ylabel('Number of Chairs (y)')
plt.title('Feasible Region and Optimal Solution')
plt.grid(True)
plt.legend()

# Create 3D plot of the profit function over the feasible region
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# Generate data for 3D plot
x_range = np.linspace(0, 30, 50)
y_range = np.linspace(0, 40, 50)
X, Y = np.meshgrid(x_range, y_range)
Z = 40*X - 0.5*X**2 + 30*Y - 0.4*Y**2

# Plot the profit surface
surf = ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8, edgecolor='none')
ax.set_xlabel('Tables (x)')
ax.set_ylabel('Chairs (y)')
ax.set_zlabel('Profit ($)')
ax.set_title('Profit Function Surface')

# Highlight the optimal point
ax.scatter([optimal_x], [optimal_y], [-objective([optimal_x, optimal_y])], 
           color='red', s=100, label='Optimal Point')

# Add a colorbar
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5, label='Profit ($)')

# Show the constraint boundaries on z=0 plane
max_z = np.max(Z)
x_wood = np.linspace(0, 30, 50)
y_wood = (300 - 10*x_wood) / 5
ax.plot(x_wood, y_wood, np.zeros_like(x_wood), 'r-', label='Wood Constraint')

x_labor = np.linspace(0, 30, 50) 
y_labor = (110 - 5*x_labor) / 3
ax.plot(x_labor, y_labor, np.zeros_like(x_labor), 'b-', label='Labor Constraint')

ax.legend()

# Show profit contours
plt.figure(figsize=(10, 6))
contour_levels = np.linspace(0, 1000, 20)
cp = plt.contour(X, Y, Z, levels=contour_levels, cmap='viridis')
plt.clabel(cp, inline=True, fontsize=8)
plt.colorbar(label='Profit ($)')

# Plot the constraints
plt.plot(x, y_wood, 'r-', label='Wood Constraint')
plt.plot(x, y_labor, 'b-', label='Labor Constraint')
plt.fill_between(x, 0, np.minimum(y_wood, y_labor), where=(x >= 0) & (np.minimum(y_wood, y_labor) >= 0), 
                 alpha=0.3, color='green', label='Feasible Region')
plt.scatter(optimal_x, optimal_y, color='red', s=100, label=f'Optimal Solution ({optimal_x:.2f}, {optimal_y:.2f})')
plt.xlim(0, 35)
plt.ylim(0, 40)
plt.xlabel('Number of Tables (x)')
plt.ylabel('Number of Chairs (y)')
plt.title('Profit Contours with Feasible Region')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

Code Explanation

Let’s break down the key components of our solution:

1. Problem Setup

We begin by importing the necessary libraries:

  • numpy for numerical operations
  • matplotlib for visualization
  • scipy.optimize for solving our nonlinear optimization problem

2. The Objective Function

Our profit function to maximize is:

$$P(x, y) = 40x - 0.5x^2 + 30y - 0.4y^2$$

But since SciPy’s minimize function minimizes rather than maximizes, we define the objective function with a negative sign:

1
2
3
def objective(vars):
    x, y = vars
    return -1 * (40*x - 0.5*x**2 + 30*y - 0.4*y**2)

3. Constraints Definition

Our constraints are:

  • Wood constraint: $10x + 5y \leq 300$
  • Labor constraint: $5x + 3y \leq 110$
  • Non-negativity: $x \geq 0, y \geq 0$

These are implemented as inequality constraints in the format SciPy expects:

1
2
3
4
5
6
constraints = [
    {'type': 'ineq', 'fun': lambda vars: 300 - 10*vars[0] - 5*vars[1]},  # Wood constraint
    {'type': 'ineq', 'fun': lambda vars: 110 - 5*vars[0] - 3*vars[1]},   # Labor constraint
    {'type': 'ineq', 'fun': lambda vars: vars[0]},  # x ≥ 0
    {'type': 'ineq', 'fun': lambda vars: vars[1]}   # y ≥ 0
]

4. Solving the Optimization Problem

We use SciPy’s minimize function with the SLSQP method (Sequential Least Squares Programming), which can handle nonlinear objective functions with constraints:

1
result = minimize(objective, initial_guess, constraints=constraints, method='SLSQP')

5. Visualizations

We create three informative visualizations:

  1. Feasible Region Plot: Shows the area where all constraints are satisfied and highlights the optimal solution
  2. 3D Profit Surface: Displays the profit function as a surface with the optimal point highlighted
  3. Profit Contours: Shows profit levels as contour lines overlaid on the feasible region

Results Analysis

Optimal solution:
Tables to produce: 12.07
Chairs to produce: 16.55
Maximum profit: $796.90

Resource utilization:
Wood used: 203.45 units out of 300 units (67.82%)
Labor used: 110.00 hours out of 110 hours (100.00%)

Let’s analyze the results of our optimization:

The optimal solution is to produce approximately 12.07 tables and 16.55 chairs, which yields a maximum profit of $796.90.

Looking at resource utilization:

  • Wood: We use almost all of our wood (203.45 units out of 300, or 67.82%)
  • Labor: We use almost all of our labor hours (110.00 hours out of 110, or 100.00%)

This tells us that both wood and labor are limiting factors in our production.
If we wanted to increase profit further, we would need to increase both resources.

Visual Insights

The visualizations provide several key insights:

  1. Feasible Region: The triangular green area shows all possible combinations of tables and chairs we could produce with our constraints.

  1. 3D Profit Surface: The curved surface illustrates how profit changes as we adjust production quantities.
    Notice that it’s not a plane (which would indicate linear profits) but has a curved shape due to the quadratic terms.

  1. Profit Contours: These lines show combinations of tables and chairs that yield the same profit.
    The optimal solution is at the point where the highest contour line touches the feasible region.

Practical Implications

This nonlinear approach to product mix offers several advantages over linear models:

  1. Market Reality: It accounts for decreasing marginal profits as production increases, which often happens due to price discounting, market saturation, or increasing marginal costs.

  2. Better Decision Making: By understanding the nonlinear nature of profits, managers can make more informed decisions about resource allocation.

  3. Sensitivity Analysis: By examining the contour plot, we can see how profit changes if we deviate slightly from the optimal solution, which helps in real-world scenarios where exact quantities might not be achievable.

Conclusion

Nonlinear programming provides a powerful approach to product mix optimization when profit functions aren’t linear.
Our Python solution demonstrates how to formulate, solve, and visualize such problems effectively.

For manufacturers, this means more realistic and profitable production planning that accounts for the complexities of real-world economics.

Optimization of Post-Disaster Infrastructure Recovery Scheduling

Today, we’ll dive into an important real-world problem: prioritizing infrastructure recovery after a disaster.
When natural disasters strike, decision-makers must determine the optimal sequence for restoring critical services with limited resources.
Let’s explore this challenge using Python optimization techniques.

The Problem: Post-Disaster Recovery Prioritization

Imagine a scenario where a major hurricane has damaged multiple infrastructure systems in a coastal city.
We need to schedule repairs to maximize the restored population service while considering:

  1. Limited repair crews/resources
  2. Different repair times for each facility
  3. Varying population impacts
  4. Dependencies between infrastructure systems

Mathematical Formulation

We can formulate this as an optimization problem.
Let’s define:

  • $i$ = infrastructure facility index
  • $p_i$ = population served by facility $i$
  • $t_i$ = time required to repair facility $i$
  • $d_{ij}$ = dependency (1 if facility $i$ depends on facility $j$, 0 otherwise)
  • $x_i$ = repair start time for facility $i$
  • $C_i$ = completion time for facility $i$ (equals $x_i + t_i$)

Our objective is to minimize the total population-hours without service:

$$\min \sum_{i} p_i \cdot C_i$$

Subject to constraints:

  • Resource constraints (limited repair crews)
  • Dependency constraints (some facilities must be repaired before others)
  • Non-negative start times

Let’s implement this using Python and solve it with PuLP, a popular linear programming library.

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
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pulp import *
import networkx as nx
from matplotlib.colors import ListedColormap

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

# Define the problem data
num_facilities = 8
facility_names = [
    "Power Plant", "Water Treatment", "Hospital", 
    "Cell Tower", "Bridge", "Shelter", 
    "Grocery Store", "Gas Station"
]

# Define repair times (days)
repair_times = {
    "Power Plant": 5,
    "Water Treatment": 3,
    "Hospital": 4,
    "Cell Tower": 2,
    "Bridge": 7,
    "Shelter": 2,
    "Grocery Store": 1,
    "Gas Station": 1
}

# Population served (thousands)
population_impact = {
    "Power Plant": 50,
    "Water Treatment": 45,
    "Hospital": 30,
    "Cell Tower": 35,
    "Bridge": 25,
    "Shelter": 20,
    "Grocery Store": 15,
    "Gas Station": 10
}

# Dependencies (which facilities depend on others)
# Format: {dependent_facility: [facilities it depends on]}
dependencies = {
    "Water Treatment": ["Power Plant"],
    "Hospital": ["Power Plant", "Water Treatment"],
    "Cell Tower": ["Power Plant"],
    "Shelter": ["Power Plant", "Water Treatment"],
    "Grocery Store": ["Power Plant", "Bridge"],
    "Gas Station": ["Power Plant", "Bridge"]
}

# Number of repair crews available
num_crews = 3

# Create a PuLP model
model = LpProblem("Disaster_Recovery_Scheduling", LpMinimize)

# Define variables
start_times = {facility: LpVariable(f"start_{i}", lowBound=0, cat='Integer') 
               for i, facility in enumerate(facility_names)}

# Maximum planning horizon (sum of all repair times as a safe upper bound)
max_horizon = sum(repair_times.values())

# Binary variables for crew assignment
# x_i_t = 1 if repair of facility i starts at time t
crew_assignments = {}
for facility in facility_names:
    for t in range(max_horizon):
        crew_assignments[(facility, t)] = LpVariable(f"assign_{facility}_{t}", cat='Binary')

# Each facility must be repaired exactly once
for facility in facility_names:
    model += lpSum(crew_assignments[(facility, t)] for t in range(max_horizon)) == 1

# Connect assignment variables with start times
for facility in facility_names:
    model += lpSum(t * crew_assignments[(facility, t)] for t in range(max_horizon)) == start_times[facility]

# Resource constraint: at most num_crews active repairs at any time
for t in range(max_horizon):
    active_repairs = []
    for facility in facility_names:
        for tau in range(max(0, t - repair_times[facility] + 1), t + 1):
            if tau < max_horizon:
                active_repairs.append(crew_assignments[(facility, tau)])
    model += lpSum(active_repairs) <= num_crews

# Dependency constraints
for facility, deps in dependencies.items():
    for dep in deps:
        # facility can only start after dep is completed
        model += start_times[facility] >= start_times[dep] + repair_times[dep]

# Calculate completion times
completion_times = {facility: start_times[facility] + repair_times[facility] for facility in facility_names}

# Objective function: minimize population impact × completion time
model += lpSum(population_impact[facility] * completion_times[facility] for facility in facility_names)

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

# Extract results
results = {
    "Facility": [],
    "Start Time": [],
    "Repair Time": [],
    "End Time": [],
    "Population Impact": []
}

for facility in facility_names:
    results["Facility"].append(facility)
    results["Start Time"].append(start_times[facility].value())
    results["Repair Time"].append(repair_times[facility])
    results["End Time"].append(start_times[facility].value() + repair_times[facility])
    results["Population Impact"].append(population_impact[facility])

results_df = pd.DataFrame(results)
results_df = results_df.sort_values("Start Time").reset_index(drop=True)
print("Optimization Status:", LpStatus[model.status])
print("\nOptimal Repair Schedule:")
print(results_df)

# Calculate the objective value
objective_value = sum(population_impact[facility] * (start_times[facility].value() + repair_times[facility]) 
                     for facility in facility_names)
print(f"\nTotal population-days without service: {objective_value:.2f} thousand person-days")

# Visualization of the schedule
plt.figure(figsize=(12, 6))
colors = plt.cm.tab10(np.linspace(0, 1, len(facility_names)))

for i, row in results_df.iterrows():
    plt.barh(row['Facility'], row['Repair Time'], left=row['Start Time'], 
             color=colors[i], alpha=0.8, edgecolor='black')
    # Add text showing population impact
    plt.text(row['Start Time'] + row['Repair Time']/2, row['Facility'], 
             f"{row['Population Impact']}k", 
             va='center', ha='center', fontweight='bold')

plt.xlabel('Days after disaster')
plt.ylabel('Infrastructure Facility')
plt.title('Optimal Recovery Schedule')
plt.grid(axis='x', linestyle='--', alpha=0.7)
plt.tight_layout()

# Create a visualization of dependencies
plt.figure(figsize=(10, 8))
G = nx.DiGraph()

# Add nodes
for facility in facility_names:
    G.add_node(facility, repair_time=repair_times[facility], 
               population=population_impact[facility],
               start_time=start_times[facility].value())

# Add edges based on dependencies
for facility, deps in dependencies.items():
    for dep in deps:
        G.add_edge(dep, facility)

# Position nodes using a hierarchical layout
pos = nx.spring_layout(G, seed=42)

# Node sizes based on population impact
node_sizes = [population_impact[facility]*20 for facility in G.nodes()]

# Node colors based on start times
start_time_values = [start_times[facility].value() for facility in G.nodes()]
node_colors = [max(start_time_values) - st for st in start_time_values]  # Reverse for better color mapping

# Create custom labels with repair times
node_labels = {facility: f"{facility}\n({repair_times[facility]}d)" for facility in G.nodes()}

# Draw the graph
nx.draw_networkx_nodes(G, pos, node_size=node_sizes, node_color=node_colors, 
                      cmap=plt.cm.viridis, alpha=0.8)
nx.draw_networkx_edges(G, pos, edge_color='gray', width=2, arrowsize=20, alpha=0.6)
nx.draw_networkx_labels(G, pos, labels=node_labels, font_size=10, font_weight='bold')

plt.title('Infrastructure Dependencies and Repair Sequence')
plt.axis('off')
plt.tight_layout()

# Create a Gantt chart with multiple crews visualization
plt.figure(figsize=(14, 8))

# Extract assignment information
assignments = {}
for facility in facility_names:
    for t in range(max_horizon):
        if crew_assignments[(facility, t)].value() > 0.5:  # Binary variable is 1
            assignments[facility] = t

# Determine which crew is assigned to each facility
crew_to_facility = {}
sorted_facilities = sorted(facility_names, key=lambda f: assignments[f])

for crew_idx in range(num_crews):
    crew_to_facility[crew_idx] = []

current_workload = [0] * num_crews

for facility in sorted_facilities:
    # Assign to crew with minimum current workload
    min_crew = current_workload.index(min(current_workload))
    crew_to_facility[min_crew].append(facility)
    current_workload[min_crew] += repair_times[facility]

# Colors for different crews
crew_colors = plt.cm.tab10(np.linspace(0, 1, num_crews))

# Plot Gantt chart by crew
current_y = 0
yticks = []
ylabels = []

for crew_idx in range(num_crews):
    crew_facilities = crew_to_facility[crew_idx]
    
    if not crew_facilities:
        continue
        
    # Add crew label
    plt.text(-1.5, current_y + len(crew_facilities)/2, f"Crew {crew_idx+1}", 
             va='center', ha='center', fontweight='bold', fontsize=12)
    
    for facility in crew_facilities:
        current_y += 1
        start = start_times[facility].value()
        duration = repair_times[facility]
        
        plt.barh(current_y, duration, left=start, height=0.6, 
                color=crew_colors[crew_idx], alpha=0.8, edgecolor='black')
        
        # Add facility name and population info
        plt.text(start + duration/2, current_y, 
                f"{facility} ({population_impact[facility]}k)", 
                va='center', ha='center', fontweight='bold')
        
        yticks.append(current_y)
        ylabels.append(facility)
    
    # Add a small gap between crews
    current_y += 0.5

plt.yticks(yticks, ylabels)
plt.xlabel('Days after disaster')
plt.title('Recovery Schedule by Crew Assignment')
plt.grid(axis='x', linestyle='--', alpha=0.7)
plt.xlim(-2, max(results_df['End Time']) + 1)
plt.tight_layout()

# Cumulative service restoration over time
def calculate_service_restored(time_point, results_df):
    completed = results_df[results_df['End Time'] <= time_point]
    return completed['Population Impact'].sum()

time_points = range(int(max(results_df['End Time'])) + 2)
service_restored = [calculate_service_restored(t, results_df) for t in time_points]
total_population = sum(population_impact.values())
percentage_restored = [100 * s / total_population for s in service_restored]

plt.figure(figsize=(12, 6))
plt.plot(time_points, service_restored, marker='o', markersize=8, linewidth=3)
plt.xlabel('Days after disaster')
plt.ylabel('Population with service restored (thousands)')
plt.title('Cumulative Service Restoration Over Time')
plt.grid(True, linestyle='--', alpha=0.7)

# Add percentage labels
for i, (x, y) in enumerate(zip(time_points, service_restored)):
    if i > 0 and service_restored[i] > service_restored[i-1]:
        plt.annotate(f"{percentage_restored[i]:.1f}%", 
                    xy=(x, y), xytext=(0, 10),
                    textcoords='offset points', ha='center',
                    fontweight='bold')

plt.tight_layout()

# Show the plots
plt.show()

Understanding the Code

Let’s break down the key components of our solution:

1. Problem Definition

We’ve created a scenario with 8 critical infrastructure facilities: Power Plant, Water Treatment, Hospital, Cell Tower, Bridge, Shelter, Grocery Store, and Gas Station. Each has:

  • A repair time (days required to complete repairs)
  • Population impact (thousands of people affected)
  • Dependencies (which facilities must be operational first)

2. Mathematical Model Implementation

The code uses PuLP, a linear programming library, to model this optimization problem:

  • Decision Variables:

    • When to start repairing each facility
    • Binary variables for crew assignments

  • Objective Function:
    Minimize total population-days without service (population × time without service)

  • Constraints:

    • Resource constraints: Limited to 3 repair crews
    • Dependency constraints: Some facilities require others to be operational first
    • Assignment constraints: Each facility must be repaired exactly once

3. Visualization

The code generates four visualizations to help understand the solution:

  1. Gantt Chart: Shows the repair schedule with start and end times
  2. Dependency Graph: Visualizes the relationships between facilities
  3. Crew Assignment Chart: Shows which crew is assigned to each facility
  4. Restoration Timeline: Tracks cumulative population service restoration

Results Analysis

The optimization model provides a clear schedule prioritizing infrastructure repairs to minimize the total impact on the population.

Key Insights

  1. Critical Path: The Power Plant is repaired first as many other facilities depend on it
  2. Resource Allocation: The 3 repair crews are allocated efficiently to avoid idle time
  3. Impact-driven: Facilities serving larger populations are generally prioritized
  4. Dependency Management: The solution respects all technical dependencies

Restoration Timeline

The cumulative service restoration graph shows how quickly we can restore services to the affected population.
With our optimized schedule:

  • 50% of service is restored within the first week
  • Nearly complete restoration occurs by day 12
  • The critical first 3-4 days focus on the most impactful infrastructure

Practical Applications

This type of optimization can be applied to various disaster recovery scenarios:

  • Hurricane/flood recovery operations
  • Earthquake response
  • Power grid restoration after major outages
  • Pandemic response for healthcare infrastructure

Decision-makers can use this approach to:

  1. Develop pre-disaster recovery plans
  2. Adapt to changing conditions during a crisis
  3. Justify resource allocation decisions
  4. Communicate expectations to affected populations

Conclusion

Infrastructure recovery scheduling is a complex optimization problem with real-world impact.
By formulating it mathematically and leveraging Python’s optimization libraries, we can create decision support tools that help communities recover faster from disasters.

The approach demonstrated here balances population needs, technical constraints, and limited resources to minimize the overall impact of infrastructure disruptions.
Similar techniques can be applied to many other resource allocation and scheduling problems in emergency management.

Resource Allocation Optimization in Multi-Agent System

A Practical Implementation

In today’s increasingly interconnected world, multi-agent systems are becoming more prevalent in various domains from economics to computer networks.
One of the most fascinating challenges in these systems is how to optimally allocate limited resources among competing agents.
Let’s dive into a concrete example and solve it step by step using Python.

The Problem: Task Allocation in a Team of Robots

Imagine we have a team of robots (our agents) that need to complete various tasks in a warehouse.
Each robot has different capabilities and efficiencies for different tasks, and we need to allocate the tasks optimally to maximize overall efficiency.

This is a classic resource allocation problem in multi-agent systems.
We’ll model this as a linear optimization problem and solve it using Python’s powerful optimization libraries.

Let’s implement this solution:

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
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linear_sum_assignment
import pandas as pd
from matplotlib.colors import LinearSegmentedColormap
from pulp import LpMaximize, LpProblem, LpVariable, lpSum, LpStatus

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

def simple_task_allocation():
    """
    A simple task allocation problem using the Hungarian algorithm
    """
    print("Example 1: Simple Task Allocation Using Hungarian Algorithm")
    print("="*60)
    
    # Define number of robots and tasks
    n_robots = 5
    n_tasks = 5
    
    # Generate random efficiency matrix (how good each robot is at each task)
    # Higher value means better efficiency
    efficiency_matrix = np.random.rand(n_robots, n_tasks)
    
    # For better visualization, scale values to be between 1 and 10
    efficiency_matrix = 1 + 9 * efficiency_matrix
    
    # Display the efficiency matrix
    print("Efficiency Matrix (Robot vs Task):")
    efficiency_df = pd.DataFrame(
        efficiency_matrix, 
        index=[f"Robot {i+1}" for i in range(n_robots)],
        columns=[f"Task {i+1}" for i in range(n_tasks)]
    )
    print(efficiency_df)
    print()
    
    # The Hungarian algorithm minimizes cost, so we need to convert efficiency to cost
    # We can do this by negating efficiency
    cost_matrix = -efficiency_matrix
    
    # Apply Hungarian algorithm
    row_ind, col_ind = linear_sum_assignment(cost_matrix)
    
    # Display the results
    total_efficiency = 0
    print("Optimal Allocation:")
    for i, j in zip(row_ind, col_ind):
        print(f"Assign Robot {i+1} to Task {j+1} (Efficiency: {efficiency_matrix[i, j]:.2f})")
        total_efficiency += efficiency_matrix[i, j]
    
    print(f"\nTotal Efficiency: {total_efficiency:.2f}")
    
    # Visualize the efficiency matrix and optimal allocation
    plt.figure(figsize=(10, 8))
    sns.heatmap(efficiency_matrix, annot=True, fmt=".2f", cmap="YlGnBu",
                xticklabels=[f"Task {i+1}" for i in range(n_tasks)],
                yticklabels=[f"Robot {i+1}" for i in range(n_robots)])
    
    # Highlight the optimal allocation
    for i, j in zip(row_ind, col_ind):
        plt.plot(j + 0.5, i + 0.5, 'ro', markersize=10)
    
    plt.title("Robot-Task Efficiency Matrix with Optimal Allocation", fontsize=14)
    plt.tight_layout()
    plt.show()

def resource_constrained_allocation():
    """
    A more complex resource allocation problem with constraints
    using linear programming
    """
    print("\nExample 2: Resource Constrained Allocation Using Linear Programming")
    print("="*60)
    
    # Define number of robots, tasks, and resource types
    n_robots = 4
    n_tasks = 6
    n_resources = 3
    
    # Generate random efficiency matrix
    efficiency_matrix = np.random.rand(n_robots, n_tasks) * 10
    
    # Generate random resource requirements for each task
    resource_requirements = np.random.randint(1, 6, size=(n_tasks, n_resources))
    
    # Generate random resource capacities for each robot
    resource_capacities = np.random.randint(5, 15, size=(n_robots, n_resources))
    
    # Display the data
    print("Task Efficiency Matrix (Robot vs Task):")
    efficiency_df = pd.DataFrame(
        efficiency_matrix, 
        index=[f"Robot {i+1}" for i in range(n_robots)],
        columns=[f"Task {i+1}" for i in range(n_tasks)]
    )
    print(efficiency_df)
    print()
    
    print("Resource Requirements for Each Task:")
    resource_req_df = pd.DataFrame(
        resource_requirements,
        index=[f"Task {i+1}" for i in range(n_tasks)],
        columns=[f"Resource {i+1}" for i in range(n_resources)]
    )
    print(resource_req_df)
    print()
    
    print("Resource Capacities for Each Robot:")
    resource_cap_df = pd.DataFrame(
        resource_capacities,
        index=[f"Robot {i+1}" for i in range(n_robots)],
        columns=[f"Resource {i+1}" for i in range(n_resources)]
    )
    print(resource_cap_df)
    print()
    
    # Create the optimization model
    model = LpProblem(name="robot-task-allocation", sense=LpMaximize)
    
    # Create decision variables
    # x[i,j] = 1 if robot i is assigned to task j, 0 otherwise
    x = {}
    for i in range(n_robots):
        for j in range(n_tasks):
            x[i, j] = LpVariable(name=f"x_{i}_{j}", cat="Binary")
    
    # Objective function: maximize total efficiency
    model += lpSum(efficiency_matrix[i, j] * x[i, j] for i in range(n_robots) for j in range(n_tasks))
    
    # Constraints:
    
    # 1. Each task can be assigned to at most one robot
    for j in range(n_tasks):
        model += lpSum(x[i, j] for i in range(n_robots)) <= 1
    
    # 2. Resource constraints for each robot
    for i in range(n_robots):
        for k in range(n_resources):
            model += lpSum(resource_requirements[j, k] * x[i, j] for j in range(n_tasks)) <= resource_capacities[i, k]
    
    # Solve the model
    model.solve()
    
    # Display results
    print(f"Status: {LpStatus[model.status]}")
    
    if LpStatus[model.status] == "Optimal":
        print("\nOptimal Allocation:")
        total_efficiency = 0
        allocation_matrix = np.zeros((n_robots, n_tasks))
        
        for i in range(n_robots):
            for j in range(n_tasks):
                if x[i, j].value() > 0.5:  # Binary variables might have small numerical errors
                    print(f"Assign Robot {i+1} to Task {j+1} (Efficiency: {efficiency_matrix[i, j]:.2f})")
                    allocation_matrix[i, j] = 1
                    total_efficiency += efficiency_matrix[i, j]
        
        print(f"\nTotal Efficiency: {total_efficiency:.2f}")
        
        # Calculate resource usage
        resource_usage = np.zeros((n_robots, n_resources))
        for i in range(n_robots):
            for j in range(n_tasks):
                if x[i, j].value() > 0.5:
                    for k in range(n_resources):
                        resource_usage[i, k] += resource_requirements[j, k]
        
        # Display resource usage
        print("\nResource Usage for Each Robot:")
        resource_usage_df = pd.DataFrame(
            resource_usage,
            index=[f"Robot {i+1}" for i in range(n_robots)],
            columns=[f"Resource {i+1}" for i in range(n_resources)]
        )
        print(resource_usage_df)
        
        # Calculate resource utilization percentage
        resource_util = resource_usage / resource_capacities * 100
        print("\nResource Utilization (%):")
        resource_util_df = pd.DataFrame(
            resource_util,
            index=[f"Robot {i+1}" for i in range(n_robots)],
            columns=[f"Resource {i+1}" for i in range(n_resources)]
        )
        print(resource_util_df)
        
        # Visualize the allocation and resource utilization
        fig, axes = plt.subplots(1, 2, figsize=(18, 8))
        
        # Plot allocation matrix
        sns.heatmap(allocation_matrix, annot=False, cmap="Blues", 
                   xticklabels=[f"Task {i+1}" for i in range(n_tasks)],
                   yticklabels=[f"Robot {i+1}" for i in range(n_robots)], ax=axes[0])
        
        # Add efficiency values where assigned
        for i in range(n_robots):
            for j in range(n_tasks):
                if allocation_matrix[i, j] > 0.5:
                    axes[0].text(j + 0.5, i + 0.5, f"{efficiency_matrix[i, j]:.1f}", 
                               ha='center', va='center', color='black', fontweight='bold')
        
        axes[0].set_title("Robot-Task Allocation", fontsize=14)
        
        # Plot resource utilization
        cmap = LinearSegmentedColormap.from_list('custom_cmap', ['white', 'lightgreen', 'green'], N=256)
        sns.heatmap(resource_util, annot=True, fmt=".1f", cmap=cmap, vmin=0, vmax=100,
                   xticklabels=[f"Resource {i+1}" for i in range(n_resources)],
                   yticklabels=[f"Robot {i+1}" for i in range(n_robots)], ax=axes[1])
        axes[1].set_title("Resource Utilization (%)", fontsize=14)
        
        plt.tight_layout()
        plt.show()
        
        # Plot task resource requirements
        plt.figure(figsize=(12, 6))
        bar_width = 0.25
        x = np.arange(n_tasks)
        
        for k in range(n_resources):
            plt.bar(x + k * bar_width, resource_requirements[:, k], 
                    width=bar_width, label=f'Resource {k+1}')
        
        plt.xlabel('Tasks')
        plt.ylabel('Resource Units Required')
        plt.title('Resource Requirements per Task')
        plt.xticks(x + bar_width, [f'Task {i+1}' for i in range(n_tasks)])
        plt.legend()
        plt.tight_layout()
        plt.show()

# Run the examples
simple_task_allocation()
resource_constrained_allocation()

Understanding the Code: A Deep Dive

Let’s break down the two approaches I’ve implemented for resource allocation in multi-agent systems:

Example 1: Simple Task Allocation with Hungarian Algorithm

The first example demonstrates a classic one-to-one matching problem where each robot can be assigned to exactly one task to maximize overall efficiency.
This is perfectly suited for the Hungarian algorithm (implemented in Python as linear_sum_assignment).

The key components of this implementation are:

  1. Efficiency Matrix: We create a matrix where each cell represents how efficiently a particular robot can perform a specific task.
    Higher values indicate better performance.

  2. Cost Matrix Conversion: Since the Hungarian algorithm minimizes cost, we negate the efficiency values to convert the maximization problem into a minimization one.

  3. Optimal Assignment: The algorithm returns the row and column indices that give us the optimal assignment.

The mathematical formulation of this problem can be represented as:

$$\max \sum_{i=1}^{n} \sum_{j=1}^{m} e_{ij} \cdot x_{ij}$$

Subject to:
$$\sum_{j=1}^{m} 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:

  • $e_{ij}$ is the efficiency of robot $i$ performing task $j$
  • $x_{ij}$ is a binary variable that equals 1 if robot $i$ is assigned to task $j$, and 0 otherwise

Example 2: Resource-Constrained Allocation with Linear Programming

The second example tackles a more complex scenario where tasks require various resources, and robots have limited resource capacities.
This requires linear programming to solve (using PuLP library).

Key components include:

  1. Multiple Resources: Each task requires different amounts of multiple resource types.

  2. Resource Constraints: Each robot has limited capacity for each resource type.

  3. Binary Variables: We use binary variables to represent whether a robot is assigned to a task.

The mathematical formulation can be expressed as:

$$\max \sum_{i=1}^{n} \sum_{j=1}^{m} e_{ij} \cdot x_{ij}$$

Subject to:
$$\sum_{i=1}^{n} x_{ij} \leq 1 \quad \forall j$$
$$\sum_{j=1}^{m} r_{jk} \cdot x_{ij} \leq c_{ik} \quad \forall i, k$$
$$x_{ij} \in {0, 1} \quad \forall i,j$$

Where:

  • $e_{ij}$ is the efficiency of robot $i$ performing task $j$
  • $x_{ij}$ is a binary variable that equals 1 if robot $i$ is assigned to task $j$
  • $r_{jk}$ is the amount of resource $k$ required for task $j$
  • $c_{ik}$ is the capacity of resource $k$ available to robot $i$

Visualizing the Results

For the first example, we created a heatmap of the efficiency matrix with red circles highlighting the optimal assignments.
This makes it easy to see which robot is matched with which task.

For the resource-constrained example, we created multiple visualizations:

  1. Allocation Matrix: Shows which robot is assigned to which task along with the efficiency values.
Example 1: Simple Task Allocation Using Hungarian Algorithm
============================================================
Efficiency Matrix (Robot vs Task):
           Task 1    Task 2    Task 3    Task 4    Task 5
Robot 1  4.370861  9.556429  7.587945  6.387926  2.404168
Robot 2  2.403951  1.522753  8.795585  6.410035  7.372653
Robot 3  1.185260  9.729189  8.491984  2.911052  2.636425
Robot 4  2.650641  3.738180  5.722808  4.887505  3.621062
Robot 5  6.506676  2.255445  3.629302  4.297257  5.104630

Optimal Allocation:
Assign Robot 1 to Task 2 (Efficiency: 9.56)
Assign Robot 2 to Task 5 (Efficiency: 7.37)
Assign Robot 3 to Task 3 (Efficiency: 8.49)
Assign Robot 4 to Task 4 (Efficiency: 4.89)
Assign Robot 5 to Task 1 (Efficiency: 6.51)

Total Efficiency: 36.82

  1. Resource Utilization: Shows how much of each robot’s resource capacity is being used (as a percentage).
Example 2: Resource Constrained Allocation Using Linear Programming
============================================================
Task Efficiency Matrix (Robot vs Task):
           Task 1    Task 2    Task 3    Task 4    Task 5    Task 6
Robot 1  7.851760  1.996738  5.142344  5.924146  0.464504  6.075449
Robot 2  1.705241  0.650516  9.488855  9.656320  8.083973  3.046138
Robot 3  0.976721  6.842330  4.401525  1.220382  4.951769  0.343885
Robot 4  9.093204  2.587800  6.625223  3.117111  5.200680  5.467103

Resource Requirements for Each Task:
        Resource 1  Resource 2  Resource 3
Task 1           5           2           2
Task 2           4           2           2
Task 3           4           4           1
Task 4           5           5           2
Task 5           5           2           1
Task 6           4           4           4

Resource Capacities for Each Robot:
         Resource 1  Resource 2  Resource 3
Robot 1          13           5          13
Robot 2          11          13          12
Robot 3           5          12          12
Robot 4           7           5          12

Status: Optimal

Optimal Allocation:
Assign Robot 1 to Task 6 (Efficiency: 6.08)
Assign Robot 2 to Task 3 (Efficiency: 9.49)
Assign Robot 2 to Task 4 (Efficiency: 9.66)
Assign Robot 3 to Task 2 (Efficiency: 6.84)
Assign Robot 4 to Task 1 (Efficiency: 9.09)

Total Efficiency: 41.16

Resource Usage for Each Robot:
         Resource 1  Resource 2  Resource 3
Robot 1         4.0         4.0         4.0
Robot 2         9.0         9.0         3.0
Robot 3         4.0         2.0         2.0
Robot 4         5.0         2.0         2.0

Resource Utilization (%):
         Resource 1  Resource 2  Resource 3
Robot 1   30.769231   80.000000   30.769231
Robot 2   81.818182   69.230769   25.000000
Robot 3   80.000000   16.666667   16.666667
Robot 4   71.428571   40.000000   16.666667

  1. Task Resource Requirements: A bar chart showing how much of each resource is required by each task.

These visualizations help us understand the trade-offs and constraints in multi-agent resource allocation problems.

Practical Applications

This type of optimization is applicable in many real-world scenarios:

  • Cloud computing resource allocation
  • Supply chain management
  • Emergency service dispatch
  • Factory production planning
  • Network bandwidth allocation

By optimizing resource allocation, we can significantly improve efficiency and reduce costs in multi-agent systems.

Conclusion

Resource allocation in multi-agent systems represents a fascinating intersection of operations research and distributed computing.
The approaches demonstrated here—Hungarian algorithm for simple matching and linear programming for constrained optimization—provide powerful tools for solving these problems.

The key takeaway is that with the right mathematical formulation and optimization techniques, we can find efficient solutions to complex resource allocation problems, even with multiple constraints and objectives.

Optimizing Hospital Operating Room Schedules with Python and Linear Programming

Optimizing Operating Room Scheduling with Python: A Practical Approach

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

🏥 The Problem: Operating Room Scheduling

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

We assume:

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

Let’s define the problem using the following notation:

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

📦 Data Setup and Model in Python

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

1
2
3
4
5
!pip install pulp
import pulp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
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
# Define surgeries
surgeries = {
    'S1': {'duration': 3, 'priority': 10},
    'S2': {'duration': 4, 'priority': 8},
    'S3': {'duration': 2, 'priority': 6},
    'S4': {'duration': 5, 'priority': 7},
    'S5': {'duration': 1, 'priority': 5}
}

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

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

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

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

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

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

# Solve the problem
prob.solve()

✅ Explanation of the Code

  1. Data Definition:

    • We define $5$ surgeries with their durations and priority scores.
    • There are $3$ available operating rooms.

  2. Model Setup:

    • We use pulp.LpProblem to define a maximization problem.
    • Binary decision variables $( x_{ij} )$ are created to indicate assignments.

  3. Objective Function:

    • We maximize the total sum of priorities for scheduled surgeries.

  4. Constraints:

    • Each surgery can be assigned to at most one OR.
    • Each OR can’t exceed its 8-hour daily limit.

  5. Solving:

    • We call .solve() to let PuLP find the optimal assignment.

📊 Visualizing the Results

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

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

df_schedule = pd.DataFrame(schedule)

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

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

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

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

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

📈 Interpretation of the Output

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

🧠 Takeaways

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

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

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

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

Optimizing Disaster Evacuation and Resource Allocation with Python

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

🆘 Optimizing Evacuation Routes and Resource Allocation During Disasters with Python

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

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

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

🧠 The Problem: Evacuation and Supply Allocation

Let’s consider a simplified city model with:

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

We’ll optimize for:

  1. Minimum total travel distance
  2. Feasible distribution of evacuees and resources

🧮 Mathematical Formulation

Let:

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

Objective:

Minimize total cost:

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

Subject to:

  1. All people in each danger zone must be evacuated:
    $$
    \sum_{j \in S} x_{ij} = p_i \quad \forall i \in D
    $$

  2. Shelter capacity must not be exceeded:
    $$
    \sum_{i \in D} x_{ij} \leq C_j \quad \forall j \in S
    $$

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

🧪 Let’s Code!

✅ Install PuLP

1
!pip install pulp

🧮 Python Code (Evacuation Optimization)

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
import pulp
import matplotlib.pyplot as plt
import networkx as nx

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

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

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

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

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

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

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

# Solve
prob.solve()

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

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

📊 Visualization

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
# Create graph for visualization
G = nx.DiGraph()

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

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

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

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

🔍 Results and Analysis

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

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

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

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

🧯 Next Steps

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

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

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

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

🛠️ Designing a Reliable Power Grid with Python

Optimizing Cost and Reliability in a Network

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

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

📘 The Problem Statement

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

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

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

🧮 Mathematical Formulation

Given:

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

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

  • $( T )$ is connected

  • The total cost $( \sum_{e \in T} c(e) )$ is minimized

  • The expected reliability is maximized:

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

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

🐍 Python Code

Let’s implement this.

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
import networkx as nx
import matplotlib.pyplot as plt

# Create a graph
G = nx.Graph()

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

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

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

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

🔍 Code Breakdown

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

📊 Results & Interpretation

Total MST Cost: 10
Expected Network Reliability: 0.7195

The algorithm produces a network that:

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

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

✅ Conclusion

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

Optimizing Investment Portfolios

Solving the Capital Allocation Problem with Python

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

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

Problem Statement

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

Mathematical Formulation

The Sharpe ratio is defined as:

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

Where:

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

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

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

And the portfolio variance is:

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

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

Let’s implement this in Python:

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
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
from matplotlib.cm import get_cmap

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

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

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

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

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

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

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

# Risk-free rate
risk_free_rate = 0.02

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

Code Explanation

Let’s break down what the code does:

  1. Setup: We first define our assets, expected returns, volatilities, and correlation matrix.
    These are the inputs to our optimization problem.

  2. Covariance Matrix: We calculate the covariance matrix from the correlation matrix and volatilities.
    This is essential for calculating portfolio risk.

  3. Portfolio Performance: We create functions to calculate the expected return and volatility of a portfolio given certain weights.

  4. Optimization: We use scipy.optimize.minimize to find the weights that maximize the Sharpe ratio (or minimize the negative Sharpe ratio).
    We set constraints (weights sum to 1) and bounds (no short selling, so weights are between $0$ and $1$).

  5. Monte Carlo Simulation: We generate $10,000$ random portfolios to visualize the efficient frontier.

  6. Performance Comparison: We compare our optimal portfolio with some alternative allocations (equal weights, stocks-heavy, bonds-heavy).

Results and Visualization

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

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

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

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

Capital Allocation for $10000:
Stocks: $1700.10
Bonds: $5480.85
Real Estate: $1113.21
Commodities: $1705.84
  1. Risk-Return Profile: The efficient frontier plot shows the relationship between risk and return for different portfolios. The optimal portfolio (marked with a red star) provides the best risk-adjusted return.

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

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

Insights and Conclusions

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

  1. Diversification Benefits: The optimal portfolio is well-diversified across all asset classes, taking advantage of the correlation structure to reduce overall risk.

  2. Risk-Return Tradeoff: There’s a clear tradeoff between risk and return, but optimization helps us find the sweet spot that maximizes risk-adjusted returns.

  3. Practical Application: With this approach, investors can make data-driven decisions about capital allocation rather than relying on intuition alone.

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

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