Optimization & Linear Programming
Optimization is the math behind making the best decisions. From resource allocation to pricing, learn to formulate and solve optimization problems.
Optimization Problem Types
Linear Programming with PuLP
from pulp import *
import numpy as np
# Manufacturing optimization
def production_optimization():
# Decision variables
x1 = LpVariable("Product_A", lowBound=0, cat='Integer')
x2 = LpVariable("Product_B", lowBound=0, cat='Integer')
x3 = LpVariable("Product_C", lowBound=0, cat='Integer')
# Problem (maximize profit)
prob = LpProblem("Manufacturing", LpMaximize)
# Objective function
prob += 40 * x1 + 30 * x2 + 50 * x3, "Total Profit"
# Constraints
prob += 2 * x1 + 1 * x2 + 3 * x3 <= 120, "Labor Hours"
prob += 4 * x1 + 2 * x2 + 1 * x3 <= 200, "Material A"
prob += 1 * x1 + 3 * x2 + 2 * x3 <= 150, "Material B"
prob += x1 + x2 + x3 <= 60, "Machine Capacity"
# Solve
prob.solve(PULP_CBC_CMD(msg=0))
print(f"Status: {LpStatus[prob.status]}")
print(f"Product A: {x1.varValue}")
print(f"Product B: {x2.varValue}")
print(f"Product C: {x3.varValue}")
print(f"Total Profit: ${value(prob.objective)}")
return {v.name: v.varValue for v in prob.variables()}
# Supply chain network design
def supply_chain_optimization():
# Facilities, warehouses, customers
facilities = ['F1', 'F2']
warehouses = ['W1', 'W2', 'W3']
customers = ['C1', 'C2', 'C3', 'C4']
# Costs: facility -> warehouse, warehouse -> customer
facility_warehouse_cost = {
('F1', 'W1'): 2, ('F1', 'W2'): 3, ('F1', 'W3'): 4,
('F2', 'W1'): 3, ('F2', 'W2'): 2, ('F2', 'W3'): 3
}
warehouse_customer_cost = {
('W1', 'C1'): 1, ('W1', 'C2'): 2, ('W1', 'C3'): 3, ('W1', 'C4'): 4,
('W2', 'C1'): 2, ('W2', 'C2'): 1, ('W2', 'C3'): 2, ('W2', 'C4'): 3,
('W3', 'C1'): 3, ('W3', 'C2'): 2, ('W3', 'C3'): 1, ('W3', 'C4'): 2
}
prob = LpProblem("SupplyChain", LpMinimize)
# Variables
fw = LpVariable.dicts("FW", facility_warehouse_cost.keys(), lowBound=0)
wc = LpVariable.dicts("WC", warehouse_customer_cost.keys(), lowBound=0)
# Objective
prob += (lpSum([facility_warehouse_cost[k] * fw[k] for k in fw]) +
lpSum([warehouse_customer_cost[k] * wc[k] for k in wc]))
# Capacity constraints
for f in facilities:
prob += lpSum([fw[k] for k in fw if k[0] == f]) <= 500
# Demand constraints
demand = {'C1': 100, 'C2': 150, 'C3': 200, 'C4': 120}
for c in customers:
prob += lpSum([wc[k] for k in wc if k[1] == c]) == demand[c]
# Flow conservation
for w in warehouses:
prob += lpSum([fw[k] for k in fw if k[1] == w]) == \
lpSum([wc[k] for k in wc if k[0] == w])
prob.solve(PULP_CBC_CMD(msg=0))
return value(prob.objective)
Convex Optimization with SciPy
from scipy.optimize import minimize, LinearConstraint, Bounds
import numpy as np
# Portfolio optimization (Markowitz)
def portfolio_optimization(returns, cov_matrix, risk_aversion=1.0):
n_assets = len(returns)
def neg_sharpe(weights):
portfolio_return = np.dot(weights, returns)
portfolio_risk = np.sqrt(np.dot(weights.T, np.dot(cov_matrix, weights)))
return -(portfolio_return - 0.02) / portfolio_risk # Risk-free rate = 2%
# Constraints
constraints = [
{'type': 'eq', 'fun': lambda w: np.sum(w) - 1} # Weights sum to 1
]
bounds = Bounds(0, 0.4) # Max 40% in any asset
# Initial guess
x0 = np.ones(n_assets) / n_assets
result = minimize(neg_sharpe, x0, method='SLSQP',
bounds=bounds, constraints=constraints)
return {
'weights': result.x,
'expected_return': np.dot(result.x, returns),
'risk': np.sqrt(np.dot(result.x.T, np.dot(cov_matrix, result.x)))
}
# Resource allocation with gradient descent
def gradient_descent_optimization():
def objective(x):
return (x[0] - 1)**2 + (x[1] - 2.5)**2
def gradient(x):
return np.array([2*(x[0] - 1), 2*(x[1] - 2.5)])
x = np.array([0.0, 0.0])
lr = 0.1
for i in range(100):
grad = gradient(x)
x = x - lr * grad
if np.linalg.norm(grad) < 1e-6:
break
return x, objective(x)
# Constrained optimization
def constrained_optimization():
def objective(x):
return x[0]**2 + x[1]**2
def objective_grad(x):
return np.array([2*x[0], 2*x[1]])
constraints = LinearConstraint(
[[1, 1], [1, -1]],
[1, -np.inf],
[np.inf, 1]
)
bounds = Bounds([0, 0], [5, 5])
result = minimize(objective, [1, 1], jac=objective_grad,
method='trust-constr',
bounds=bounds, constraints=constraints)
return result.x, result.fun
Integer Programming for Scheduling
from pulp import *
import pandas as pd
def employee_scheduling():
# Requirements per shift
shifts = ['Morning', 'Afternoon', 'Night']
days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
requirements = {
'Morning': [5, 5, 5, 5, 5, 3, 3],
'Afternoon': [4, 4, 4, 4, 4, 2, 2],
'Night': [2, 2, 2, 2, 2, 1, 1]
}
n_employees = 15
prob = LpProblem("Scheduling", LpMinimize)
# Binary: employee i works shift j on day k
x = LpVariable.dicts("x",
((i, j, k) for i in range(n_employees) for j in shifts for k in days),
cat='Binary')
# Minimize total shifts
prob += lpSum(x.values())
# Meet requirements
for j in shifts:
for k_idx, k in enumerate(days):
prob += lpSum(x[i, j, k] for i in range(n_employees)) >= requirements[j][k_idx]
# Each employee works max 5 days
for i in range(n_employees):
prob += lpSum(x[i, j, k] for j in shifts for k in days) <= 5
# At most one shift per day
for i in range(n_employees):
for k in days:
prob += lpSum(x[i, j, k] for j in shifts) <= 1
prob.solve(PULP_CBC_CMD(msg=0))
schedule = []
for i in range(n_employees):
for j in shifts:
for k in days:
if value(x[i, j, k]) == 1:
schedule.append({'Employee': i, 'Shift': j, 'Day': k})
return pd.DataFrame(schedule)
# Vehicle routing
def vehicle_routing():
n_locations = 10
n_vehicles = 3
np.random.seed(42)
distances = np.random.randint(1, 50, (n_locations, n_locations))
np.fill_diagonal(distances, 0)
prob = LpProblem("VRP", LpMinimize)
x = LpVariable.dicts("x",
((i, j, k) for i in range(n_locations)
for j in range(n_locations)
for k in range(n_vehicles)),
cat='Binary')
u = LpVariable.dicts("u", range(1, n_locations), lowBound=1, upBound=n_locations-1)
# Minimize total distance
prob += lpSum(distances[i,j] * x[i,j,k]
for i in range(n_locations)
for j in range(n_locations)
for k in range(n_vehicles))
# Each location visited exactly once
for j in range(1, n_locations):
prob += lpSum(x[i,j,k] for i in range(n_locations)
for k in range(n_vehicles)) == 1
# Enter and leave each location
for k in range(n_vehicles):
for i in range(n_locations):
prob += lpSum(x[i,j,k] for j in range(n_locations)) == \
lpSum(x[j,i,k] for j in range(n_locations))
# Subtour elimination
for k in range(n_vehicles):
for i in range(1, n_locations):
for j in range(1, n_locations):
if i != j:
prob += u[i] - u[j] + n_locations * x[i,j,k] <= n_locations - 1
prob.solve(PULP_CBC_CMD(msg=0))
return value(prob.objective)
Best Practices
- Formulate carefully β small modeling errors lead to wrong solutions
- Start with LP relaxations to understand the problem structure
- Use warm starting for iterative optimization
- Validate solutions against business constraints
- Consider decomposition for very large problems