Contents  |  ‹ Previous  |  Next ›  |  Download EPUB  |  PDF

Chapter 12
Software - Python

Learning Outcomes

12.1 Why Python for Optimization?

Python has become the dominant programming language for data science, machine learning, and increasingly, operations research and optimization. But why Python? After all, Python is an interpreted language, meaning it runs slower than compiled languages like C++ or Java. The answer lies in Python’s role as a glue language and its incredible ecosystem.

12.1.1 The Power of Python

When you solve an optimization problem in Python using a library like PuLP, here’s what actually happens:

1.
You write your model in Python’s readable, intuitive syntax
2.
Python translates your model into a standard format (like LP or MPS files)
3.
Python calls a highly optimized solver written in C/C++ (like CBC, CPLEX, or Gurobi)
4.
The solver does the heavy computational work at near-optimal speed
5.
Python receives the results and lets you analyze them with powerful data tools

This means you get the best of both worlds: the ease and flexibility of Python for modeling and analysis, combined with the raw computational power of industrial-strength solvers.

Info

Python’s Superpowers for Optimization:

Compare writing a constraint in Python:

prob += lpSum(cost[i,j] * x[i,j] for i in plants for j in markets) <= budget

To writing the same thing in a lower-level language – Python wins hands down for clarity and development speed.

12.1.2 Jupyter Notebooks: A Data Analyst’s Best Friend

For this book, we’ll primarily use Jupyter Notebooks. If you haven’t encountered them before, Jupyter notebooks are interactive documents that combine:

This interactive style is perfect for optimization work where you want to:

Warning

When Notebooks Aren’t Ideal:

Jupyter notebooks excel at exploration and analysis, but they’re not ideal for:

For serious software development, you’d move your polished code into .py files and use an IDE like VS Code or PyCharm. But for learning, prototyping, and one-off analyses, notebooks are hard to beat.

12.2 Getting Started

12.2.1 Options for Running Python

You have several options for running Python and Jupyter notebooks:

Google Colab (Recommended for Beginners) Google Colab (https://colab.research.google.com) is a free, cloud-based Jupyter environment. Benefits:

To get started:

1.
Go to https://colab.research.google.com
2.
Sign in with your Google account
3.
Click “New Notebook” to create a blank notebook
4.
Start coding!

Anaconda (Recommended for Local Installation) Anaconda (https://www.anaconda.com) is a Python distribution designed for data science. It includes:

After installing Anaconda, launch “Jupyter Notebook” from the Anaconda Navigator or by typing jupyter notebook in your terminal.

VS Code Visual Studio Code (https://code.visualstudio.com) is a popular code editor that supports Jupyter notebooks through extensions. It’s a good choice if you want a more traditional development environment but still want notebook capabilities.

12.2.2 Installing Packages

Most packages we need come pre-installed in Colab or Anaconda. The main package you’ll need to install is PuLP for optimization modeling:

# In a Jupyter cell, use ! to run shell commands 
!pip install pulp

Important: After installing a new package, you need to restart your notebook kernel for the changes to take effect. In Jupyter, go to Kernel → Restart. In Colab, go to Runtime → Restart runtime.

Once installed, import the packages you need at the top of your notebook:

# Standard imports for optimization work 
import numpy as np # Numerical computing 
import pandas as pd # Data manipulation 
import matplotlib.pyplot as plt # Visualization 
from pulp import * # Optimization modeling

Resources

Interactive Practice Notebooks

Reference Materials

Try it out visually!

One LP, Four Modeling Languages: the same problem in PuLP, Pyomo, gurobipy, and AMPL with a code stepper.

12.3 Python Basics

Before diving into optimization-specific Python, let’s cover the fundamental building blocks of the language. If you’re completely new to programming, this section will get you started. If you have experience with other languages, you’ll see how Python’s syntax differs.

Resources

A companion Jupyter notebook intro_to_python_extended.ipynb is available in the book repository’s code/python/ folder. We recommend working through this notebook interactively as you read this section.

12.3.1 Variables and Data Types

Python is dynamically typed, meaning you don’t need to declare variable types explicitly; Python figures it out from context. The basic data types are:

# Variable assignment - no type declarations needed 
x = 10 # Integer 
y = 3.14 # Float 
name = "Alice" # String 
is_happy = True # Boolean 
 
print(x, y, name, is_happy) 
# Output: 10 3.14 Alice True

12.3.2 Arithmetic Operations

Python supports all standard arithmetic operations:

a = 15 
b = 4 
 
print("Addition:", a + b) # 19 
print("Subtraction:", a - b) # 11 
print("Multiplication:", a * b) # 60 
print("Division:", a / b) # 3.75 
print("Floor division:", a // b) # 3 (rounds down) 
print("Modulus:", a % b) # 3 (remainder) 
print("Exponentiation:", a ** b) # 50625 (15^4)

Info

Division Tip: In Python 3, / always returns a float, even for integers. Use // for integer division that rounds down to the nearest whole number.

12.3.3 Collections: Lists, Tuples, Dictionaries, and Sets

Python provides several built-in data structures for organizing data:

# Lists - mutable, ordered sequences 
my_list = [1, 2, 3, 4] 
 
# Tuples - immutable sequences (can't be changed) 
my_tuple = (10, 20, 30) 
 
# Dictionaries - key-value pairs 
my_dict = {"apple": 2, "banana": 3} 
 
# Sets - unique elements only (duplicates removed) 
my_set = {1, 2, 2, 3, 4} # becomes {1, 2, 3, 4} 
 
print("List:", my_list) 
print("Tuple:", my_tuple) 
print("Dictionary:", my_dict) 
print("Set:", my_set)

You can modify lists and dictionaries after creation:

# List operations 
my_list.append(5) # Add element: [1, 2, 3, 4, 5] 
my_list.remove(2) # Remove element: [1, 3, 4, 5] 
my_list[0] = 100 # Change element: [100, 3, 4, 5] 
 
# Dictionary operations 
my_dict["cherry"] = 5 # Add new key-value pair 
my_dict["apple"] = 10 # Update existing value 
print(my_dict["banana"]) # Access value by key: 3

12.3.4 Indexing and Slicing

Python uses zero-based indexing—the first element is at index 0, not 1. You can also use negative indices to count from the end.

colors = ["red", "green", "blue", "yellow", "purple"] 
 
# Single element access 
print(colors[0]) # "red" (first element) 
print(colors[-1]) # "purple" (last element) 
print(colors[-2]) # "yellow" (second to last) 
 
# Slicing: [start:end] - end is exclusive! 
print(colors[1:3]) # ["green", "blue"] 
print(colors[:3]) # ["red", "green", "blue"] (from beginning) 
print(colors[2:]) # ["blue", "yellow", "purple"] (to end) 
 
# With step: [start:end:step] 
print(colors[::2]) # ["red", "blue", "purple"] (every 2nd)

Warning

Common Mistake: Remember that slicing with [a:b] includes index a but excludes index b. So colors[1:3] gives you elements at indices 1 and 2, not 1, 2, and 3.

12.3.5 Control Flow: Conditionals

Python uses indentation (not braces) to define code blocks. The if, elif, and else keywords handle conditional logic:

num = 7 
 
if num % 2 == 0: 
   print(f"{num} is even.") 
else: 
   print(f"{num} is odd.") 
# Output: 7 is odd. 
 
# Multiple conditions with elif 
score = 85 
if score >= 90: 
   grade = "A" 
elif score >= 80: 
   grade = "B" 
elif score >= 70: 
   grade = "C" 
else: 
   grade = "F" 
print(f"Score {score} earns grade {grade}") # Grade B

Info

f-strings: The f"..." syntax creates formatted strings. Variables inside curly braces {var} are automatically converted to text. This is much cleaner than string concatenation.

12.3.6 Control Flow: Loops

For Loops For loops iterate over sequences (lists, strings, ranges, etc.):

# Loop over a list 
fruits = ["apple", "banana", "cherry"] 
for fruit in fruits: 
   print(f"I like {fruit}") 
 
# Loop over a range of numbers 
for i in range(5): 
   print(i) # Prints 0, 1, 2, 3, 4 
 
# Loop over dictionary items 
prices = {"apple": 1.50, "banana": 0.75, "cherry": 2.00} 
for fruit, price in prices.items(): 
   print(f"{fruit} costs ${price:.2f}")

While Loops While loops continue as long as a condition is true:

count = 0 
while count < 5: 
   print(count) 
   count += 1 # Don't forget to update the condition! 
# Prints 0, 1, 2, 3, 4

12.3.7 Functions

Functions are reusable blocks of code defined with the def keyword:

def square(x): 
   """Return the square of x.""" 
   return x * x 
 
val = 5 
print(f"The square of {val} is {square(val)}.") 
# Output: The square of 5 is 25.

Functions can have multiple parameters and default values:

def greet(name, greeting="Hello"): 
   """Greet someone with a customizable greeting.""" 
   return f"{greeting}, {name}!" 
 
print(greet("Alice")) # "Hello, Alice!" 
print(greet("Bob", "Hi there")) # "Hi there, Bob!"

Info

Docstrings: The triple-quoted string immediately after def is called a docstring. It documents what the function does and appears when you use help(function_name).

12.4 Python Essentials for Optimization

Now that we’ve covered the basics, let’s see how these Python features are used specifically in optimization modeling. The following patterns will appear repeatedly in your LP and IP work.

12.4.1 Lists: Ordered Collections

Lists are Python’s workhorse data structure – ordered collections that can hold any type of data. In optimization, we use lists to represent:

# Define sets for a transportation problem 
plants = ['Chicago', 'Denver', 'Atlanta'] 
markets = ['NYC', 'LA', 'Houston', 'Miami'] 
 
# Numerical data as lists 
supply = [100, 150, 200] 
demand = [80, 120, 90, 110]

Lists are zero-indexed, meaning the first element is at position 0:

plants[0] # 'Chicago' (first element) 
plants[1] # 'Denver' (second element) 
plants[-1] # 'Atlanta' (last element) 
len(plants) # 3 (number of elements)

Common operations you’ll use constantly:

sum(supply) # 450 (total supply) 
max(demand) # 120 (largest demand) 
min(demand) # 80 (smallest demand) 
plants.append('Boston') # Add to end of list 
'Denver' in plants # True (membership test)

12.4.2 Dictionaries: Key-Value Mappings

Dictionaries map keys to values. They’re the natural way to represent indexed parameters in optimization – think of them as lookup tables.

# Supply indexed by plant name 
supply = { 
   'Chicago': 100, 
   'Denver': 150, 
   'Atlanta': 200 
} 
 
# Access values using keys 
supply['Chicago'] # 100 
supply['Denver'] # 150

For parameters indexed by multiple sets (like transportation costs indexed by origin AND destination), use tuple keys:

# Transportation cost from plant to market 
cost = { 
   ('Chicago', 'NYC'): 3.5, 
   ('Chicago', 'LA'): 4.8, 
   ('Chicago', 'Houston'): 2.1, 
   ('Denver', 'NYC'): 4.2, 
   ('Denver', 'LA'): 2.1, 
   ('Denver', 'Houston'): 3.0, 
   # ... and so on 
} 
 
# Access using tuple key 
cost[('Chicago', 'NYC')] # 3.5 
cost[('Denver', 'LA')] # 2.1

This pattern – dictionaries with tuple keys – is fundamental to optimization modeling in Python. It directly mirrors the mathematical notation cij for a cost from origin i to destination j.

Useful dictionary operations:

supply.keys() # All keys: dict_keys(['Chicago', 'Denver', 'Atlanta']) 
supply.values() # All values: dict_values([100, 150, 200]) 
supply.items() # Key-value pairs (for iteration) 
sum(supply.values()) # 450 (total supply)

12.4.3 For Loops: Iteration

Loops let you repeat operations over collections. In optimization, you’ll use them to:

# Basic loop over a list 
for plant in plants: 
   print(f"Processing plant: {plant}")

Output: 

Processing plant: Chicago
Processing plant: Denver
Processing plant: Atlanta

When you need both the index and value, use enumerate:

for i, plant in enumerate(plants): 
   print(f"Plant {i}: {plant}")

For dictionaries, iterate over key-value pairs:

for plant, capacity in supply.items(): 
   print(f"{plant} can supply {capacity} units")

Nested loops create all combinations – essential for generating variables and constraints over multiple index sets:

# All plant-to-market routes 
for plant in plants: 
   for market in markets: 
      print(f"Route: {plant} -> {market}")

This creates |plants|×|markets| = 3 × 4 = 12 route combinations.

12.4.4 List Comprehensions: Concise List Creation

List comprehensions are a compact way to create lists. They’re Python’s way of saying “give me a list of X for each Y in Z.”

# Create list of all routes (same as nested loop above, but one line) 
routes = [(p, m) for p in plants for m in markets]

This produces: 

[(’Chicago’, ’NYC’), (’Chicago’, ’LA’), (’Chicago’, ’Houston’),
 (’Chicago’, ’Miami’), (’Denver’, ’NYC’), ... ]

Add conditions with if:

# Only plants with capacity > 120 
large_plants = [p for p, cap in supply.items() if cap > 120] 
# Result: ['Denver', 'Atlanta']

You’ll see comprehensions used extensively in PuLP for creating sums:

# Total cost (will make sense after PuLP section) 
total_cost = lpSum([cost[p,m] * x[p,m] for p in plants for m in markets])

12.4.5 Functions: Reusable Code

Functions package code for reuse. Define them with def:

def calculate_total_cost(flows, costs): 
   """ 
   Calculate total transportation cost. 
 
   Args: 
      flows: dict mapping (plant, market) to flow amount 
      costs: dict mapping (plant, market) to unit cost 
 
   Returns: 
      Total cost as a float 
   """ 
   total = 0 
   for route, amount in flows.items(): 
      total += amount * costs[route] 
   return total

The triple-quoted string is a docstring – documentation that explains what the function does. Good practice!

A simpler version using sum and a comprehension:

def calculate_total_cost(flows, costs): 
   """Calculate total transportation cost.""" 
   return sum(flows[r] * costs[r] for r in flows)

12.5 NumPy: Numerical Computing

NumPy (Numerical Python) provides efficient array operations. While PuLP doesn’t require NumPy, it’s invaluable for:

import numpy as np

12.5.1 Creating Arrays

NumPy arrays are like lists, but optimized for numerical operations:

# From a list 
costs = np.array([3, 5, 2, 4]) 
 
# 2D array (matrix) 
A = np.array([ 
   [2, 1, 3], 
   [1, 2, 1], 
   [3, 1, 2] 
]) 
 
# Useful shortcuts 
zeros = np.zeros(5) # [0, 0, 0, 0, 0] 
ones = np.ones(3) # [1, 1, 1] 
range_arr = np.arange(10) # [0, 1, 2, ..., 9]

12.5.2 Array Operations

NumPy shines at element-wise and aggregate operations:

costs = np.array([3, 5, 2, 4]) 
 
# Aggregate operations 
costs.sum() # 14 
costs.mean() # 3.5 
costs.min() # 2 
costs.max() # 5 
costs.std() # Standard deviation 
 
# Element-wise operations 
costs * 2 # array([6, 10, 4, 8]) 
costs + 1 # array([4, 6, 3, 5]) 
costs > 3 # array([False, True, False, True])

12.5.3 Linear Algebra

For matrix operations (useful in understanding LP theory):

A = np.array([[2, 1], [1, 3]]) 
b = np.array([4, 5]) 
c = np.array([3, 2]) 
x = np.array([1, 2]) 
 
# Matrix-vector product: Ax 
np.dot(A, x) # array([4, 7]) 
 
# Dot product: c'x 
np.dot(c, x) # 7 
 
# Matrix properties 
A.shape # (2, 2) 
A.T # Transpose 
np.linalg.det(A) # Determinant 
np.linalg.inv(A) # Inverse (if exists)

12.5.4 Visualizing Array Operations

NumPy array operations can be confusing at first, especially when working with slicing, stacking, and axis-specific operations. The following visual guide demonstrates key concepts that will help you understand how NumPy manipulates arrays, skills you’ll use when preparing data for optimization models.

12.6 Pandas: Data Management

Pandas is the go-to library for working with tabular data. It excels at:

import pandas as pd

12.6.1 DataFrames: Tables of Data

A DataFrame is like a spreadsheet or database table:

# Create from a dictionary 
plants_df = pd.DataFrame({ 
   'Plant': ['Chicago', 'Denver', 'Atlanta'], 
   'Capacity': [100, 150, 200], 
   'FixedCost': [10000, 12000, 8000] 
}) 
 
print(plants_df)

Output: 

     Plant  Capacity  FixedCost
0  Chicago       100      10000
1   Denver       150      12000
2  Atlanta       200       8000

12.6.2 Reading Data from Files

Real optimization problems get data from files, not hardcoded values:

# Read CSV file 
df = pd.read_csv('transportation_data.csv') 
 
# Read Excel file 
df = pd.read_excel('problem_data.xlsx', sheet_name='Costs') 
 
# Preview the data 
df.head() # First 5 rows 
df.info() # Column types and counts 
df.describe() # Summary statistics

12.6.3 Accessing and Filtering Data

# Access a column 
plants_df['Capacity'] 
 
# Filter rows 
large_plants = plants_df[plants_df['Capacity'] > 120] 
 
# Multiple conditions 
selected = plants_df[(plants_df['Capacity'] > 100) & 
               (plants_df['FixedCost'] < 11000)]

12.6.4 From DataFrame to Optimization Data

A critical skill is converting DataFrames to the dictionaries PuLP expects:

# DataFrame to dictionary 
supply_dict = dict(zip(plants_df['Plant'], plants_df['Capacity'])) 
# Result: {'Chicago': 100, 'Denver': 150, 'Atlanta': 200} 
 
# For cost matrices, you might have data like: 
# costs_df with columns: Origin, Destination, Cost 
cost_dict = {(row['Origin'], row['Destination']): row['Cost'] 
         for _, row in costs_df.iterrows()}

12.6.5 Analyzing Results

After solving, put results back into DataFrames for analysis:

# Create results DataFrame 
results = pd.DataFrame([ 
   {'Route': f"{p}->{m}", 'Flow': x[p,m].varValue, 'Cost': cost[p,m]} 
   for p in plants for m in markets 
   if x[p,m].varValue > 0 
]) 
 
# Analyze 
print(results) 
print(f"Total flow: {results['Flow'].sum()}") 
print(f"Average cost: {results['Cost'].mean():.2f}")

12.7 Matplotlib: Visualization

Visualization helps you understand your data before optimization and communicate results after. Matplotlib is Python’s foundational plotting library.

import matplotlib.pyplot as plt

12.7.1 Bar Charts: Comparing Categories

Perfect for showing supply, demand, or solution values by category:

plants = ['Chicago', 'Denver', 'Atlanta'] 
supply = [100, 150, 200] 
 
plt.figure(figsize=(8, 5)) 
plt.bar(plants, supply, color='steelblue', edgecolor='black') 
plt.xlabel('Plant') 
plt.ylabel('Supply Capacity') 
plt.title('Supply Capacity by Plant') 
plt.grid(axis='y', alpha=0.3) 
plt.show()

12.7.2 Grouped Bar Charts: Comparing Multiple Series

import numpy as np 
 
plants = ['Chicago', 'Denver', 'Atlanta'] 
supply = [100, 150, 200] 
used = [95, 150, 170] # After optimization 
 
x = np.arange(len(plants)) 
width = 0.35 
 
fig, ax = plt.subplots(figsize=(8, 5)) 
ax.bar(x - width/2, supply, width, label='Capacity', color='steelblue') 
ax.bar(x + width/2, used, width, label='Used', color='orange') 
 
ax.set_xlabel('Plant') 
ax.set_ylabel('Units') 
ax.set_title('Capacity vs. Utilization') 
ax.set_xticks(x) 
ax.set_xticklabels(plants) 
ax.legend() 
plt.show()

12.7.3 Heatmaps: Visualizing Matrices

Great for cost matrices or flow solutions:

# Cost matrix as 2D array 
costs = np.array([ 
   [3.5, 4.8, 2.1, 3.9], 
   [4.2, 2.1, 3.0, 4.5], 
   [2.8, 5.1, 2.9, 1.7] 
]) 
 
plants = ['Chicago', 'Denver', 'Atlanta'] 
markets = ['NYC', 'LA', 'Houston', 'Miami'] 
 
plt.figure(figsize=(8, 6)) 
plt.imshow(costs, cmap='YlOrRd', aspect='auto') 
plt.colorbar(label='Cost per Unit') 
plt.xticks(range(len(markets)), markets) 
plt.yticks(range(len(plants)), plants) 
plt.xlabel('Destination') 
plt.ylabel('Origin') 
plt.title('Transportation Cost Matrix') 
 
# Add text annotations 
for i in range(len(plants)): 
   for j in range(len(markets)): 
      plt.text(j, i, f'{costs[i,j]:.1f}', ha='center', va='center') 
 
plt.tight_layout() 
plt.show()

12.7.4 Plot Customization Reference

Matplotlib offers extensive customization options for creating publication-quality figures. The following guide provides a reference for colors, line styles, markers, and other formatting options that you’ll find useful when visualizing optimization results.

12.8 NetworkX: Graph Analysis

NetworkX is Python’s premier library for working with graphs and networks. It’s essential for the discrete algorithms we’ll cover in Part II, and useful for visualizing network-based optimization problems like transportation and facility location.

import networkx as nx

12.8.1 Why Graphs Matter for Optimization

Many optimization problems have an underlying network structure:

NetworkX lets you represent, analyze, and visualize these structures.

12.8.2 Creating Graphs

# Create a directed graph (edges have direction) 
G = nx.DiGraph() 
 
# Add nodes (can include attributes) 
G.add_node('Chicago', node_type='plant', supply=100) 
G.add_node('Denver', node_type='plant', supply=150) 
G.add_node('NYC', node_type='market', demand=80) 
G.add_node('LA', node_type='market', demand=120) 
 
# Add edges with attributes 
G.add_edge('Chicago', 'NYC', cost=3.5, capacity=100) 
G.add_edge('Chicago', 'LA', cost=4.8, capacity=80) 
G.add_edge('Denver', 'NYC', cost=4.2, capacity=90) 
G.add_edge('Denver', 'LA', cost=2.1, capacity=150)

12.8.3 Accessing Graph Data

# Basic properties 
G.number_of_nodes() # 4 
G.number_of_edges() # 4 
list(G.nodes()) # ['Chicago', 'Denver', 'NYC', 'LA'] 
list(G.edges()) # [('Chicago', 'NYC'), ...] 
 
# Node attributes 
G.nodes['Chicago']['supply'] # 100 
 
# Edge attributes 
G.edges['Chicago', 'NYC']['cost'] # 3.5 
 
# Neighbors (outgoing edges for DiGraph) 
list(G.successors('Chicago')) # ['NYC', 'LA']

12.8.4 Graph Algorithms

NetworkX includes many useful algorithms:

# Shortest path (by cost) 
path = nx.shortest_path(G, 'Chicago', 'LA', weight='cost') 
length = nx.shortest_path_length(G, 'Chicago', 'LA', weight='cost') 
print(f"Shortest path: {path}, cost: {length}") 
 
# All shortest paths from one source 
all_paths = nx.single_source_shortest_path(G, 'Chicago') 
 
# Check connectivity 
nx.is_weakly_connected(G) # True if underlying undirected graph is connected

12.8.5 Visualization

Seeing your network helps understand problem structure:

plt.figure(figsize=(10, 6)) 
 
# Layout algorithm positions nodes 
pos = nx.spring_layout(G, seed=42) 
 
# Draw nodes 
nx.draw_networkx_nodes(G, pos, node_color='lightblue', 
                 node_size=1500, alpha=0.9) 
 
# Draw edges 
nx.draw_networkx_edges(G, pos, edge_color='gray', 
                 arrows=True, arrowsize=20) 
 
# Draw labels 
nx.draw_networkx_labels(G, pos, font_size=10) 
 
# Edge labels (costs) 
edge_labels = nx.get_edge_attributes(G, 'cost') 
nx.draw_networkx_edge_labels(G, pos, edge_labels, font_size=8) 
 
plt.title('Transportation Network') 
plt.axis('off') 
plt.tight_layout() 
plt.show()

Info

NetworkX is particularly important for Part II of this book, where we cover graph algorithms including shortest paths, minimum spanning trees, and network flows.

12.9 GeoPandas: Geographic Visualization

Many optimization problems have a geographic component – facility locations, delivery routes, service territories. GeoPandas extends Pandas with geographic capabilities, making your models and results more tangible and easier to communicate.

import geopandas as gpd 
from shapely.geometry import Point

12.9.1 Why Geographic Visualization?

A transportation problem becomes much more intuitive when you can see the facilities and customers on a map. Benefits include:

12.9.2 Creating Geographic Data

# Facility data with coordinates 
facilities = pd.DataFrame({ 
   'Name': ['Chicago Plant', 'Denver Plant', 'Atlanta Plant'], 
   'Lat': [41.88, 39.74, 33.75], 
   'Lon': [-87.63, -104.99, -84.39], 
   'Capacity': [100, 150, 200] 
}) 
 
# Convert to GeoDataFrame 
geometry = [Point(lon, lat) for lon, lat in zip(facilities['Lon'], facilities['Lat'])] 
facilities_gdf = gpd.GeoDataFrame(facilities, geometry=geometry, crs='EPSG:4326')

12.9.3 Map Visualization

fig, ax = plt.subplots(figsize=(12, 8)) 
 
# Plot facilities (sized by capacity) 
facilities_gdf.plot(ax=ax, color='red', markersize=facilities_gdf['Capacity'], 
               alpha=0.6, edgecolor='black') 
 
# Add labels 
for idx, row in facilities_gdf.iterrows(): 
   ax.annotate(row['Name'], 
            xy=(row.geometry.x, row.geometry.y), 
            xytext=(5, 5), textcoords='offset points', 
            fontsize=9, fontweight='bold') 
 
ax.set_xlabel('Longitude') 
ax.set_ylabel('Latitude') 
ax.set_title('Facility Locations (size = capacity)') 
ax.grid(True, alpha=0.3) 
plt.show()

Info

GeoPandas is optional for this book, but it transforms abstract optimization problems into concrete, visual stories. It’s especially valuable when presenting results to stakeholders who may not be familiar with mathematical notation.

12.10 PuLP: Optimization Modeling

Now we arrive at the main event: PuLP, a Python library for formulating and solving linear programs (LP) and integer programs (IP). PuLP provides a clean, Pythonic interface for expressing optimization models that reads almost like mathematical notation.

12.10.1 Installing PuLP

!pip install pulp

After installing, restart your kernel, then import:

from pulp import *

PuLP comes bundled with the CBC solver (COIN-OR Branch and Cut), a capable open-source solver. For larger problems, you can connect PuLP to commercial solvers like Gurobi or CPLEX.

12.10.2 The PuLP Modeling Pattern

Every PuLP model follows six steps:

1.
Create a problem object
2.
Define decision variables
3.
Set the objective function
4.
Add constraints
5.
Solve the problem
6.
Extract and analyze results

Let’s work through a complete example.

12.10.3 Example: Product Mix Problem

A company makes two products. Each unit of Product 1 yields $3 profit; each unit of Product 2 yields $2 profit. Production is limited by three resources:

Resource Product 1 Product 2 Available
Material 10 5 300
Labor 4 4 160
Machine 2 6 180
Table 12.1: Resource usage per unit and availability for the two products.

Model:

max ⁡ 3X1 + 2X2  s.t.  10X1 + 5X2 ≤ 300 4X1 + 4X2 ≤ 160 2X1 + 6X2 ≤ 180 X1,X2 ≥ 0

Step 1: Create the Problem

# Create a maximization problem 
prob = LpProblem("Product_Mix", LpMaximize)
The first argument is a name (used in output files), the second specifies the sense: LpMaximize or LpMinimize.

Step 2: Define Decision Variables

# Create non-negative continuous variables 
x1 = LpVariable("X1", lowBound=0) 
x2 = LpVariable("X2", lowBound=0)
The lowBound=0 enforces X1,X2 ≥ 0. By default, variables are continuous.

Step 3: Set the Objective Function

# Add objective function 
prob += 3*x1 + 2*x2, "Total_Profit"
The += operator adds to the problem. The second argument is an optional name.

Step 4: Add Constraints

# Add constraints 
prob += 10*x1 + 5*x2 <= 300, "Material" 
prob += 4*x1 + 4*x2 <= 160, "Labor" 
prob += 2*x1 + 6*x2 <= 180, "Machine"
Each constraint uses <=, >=, or ==. The names help identify constraints in output.

Warning

Common Mistake: If you forget the comparison operator (<=, >=, ==), you’ll overwrite the objective function instead of adding a constraint! Recent versions of PuLP print a UserWarning when this happens, but the model still solves, so the mistake is easy to miss.

Wrong:

prob += 10*x1 + 5*x2 # Overwrites objective!

Correct:

prob += 10*x1 + 5*x2 <= 300 # Adds constraint

Step 5: Solve

# Solve the problem 
prob.solve() 
 
# Check status 
print(f"Status: {LpStatus[prob.status]}")
Status will be one of: Optimal, Infeasible, Unbounded, Not Solved.

Step 6: Extract Results

# Optimal objective value 
print(f"Maximum Profit: ${value(prob.objective)}") 
 
# Optimal variable values 
print(f"X1 = {x1.varValue}") 
print(f"X2 = {x2.varValue}") 
 
# Loop through all variables 
for v in prob.variables(): 
   print(f"{v.name} = {v.varValue}")
Output: 
Status: Optimal
Maximum Profit: $100.0
X1 = 20.0
X2 = 20.0

12.10.4 Variable Types

PuLP supports three variable types:

# Continuous (default) - for LP 
x = LpVariable("x", lowBound=0) 
 
# Integer - for IP 
y = LpVariable("y", lowBound=0, cat='Integer') 
 
# Binary (0 or 1) - for yes/no decisions 
z = LpVariable("z", cat='Binary')

12.10.5 Creating Multiple Variables

For problems with many variables, create them in bulk using LpVariable.dicts:

plants = ['Chicago', 'Denver', 'Atlanta'] 
markets = ['NYC', 'LA', 'Houston'] 
 
# Create all route combinations 
routes = [(p, m) for p in plants for m in markets] 
 
# Create a dictionary of variables 
x = LpVariable.dicts("Ship", routes, lowBound=0, cat='Continuous')

Now x[(’Chicago’, ’NYC’)] is a variable representing shipments from Chicago to NYC.

12.10.6 The lpSum Function

For sums over index sets, use lpSum (more efficient than Python’s sum):

# Objective: minimize total cost 
prob += lpSum([cost[p,m] * x[p,m] for p in plants for m in markets]) 
 
# Supply constraint for each plant 
for p in plants: 
   prob += lpSum([x[p,m] for m in markets]) <= supply[p], f"Supply_{p}" 
 
# Demand constraint for each market 
for m in markets: 
   prob += lpSum([x[p,m] for p in plants]) >= demand[m], f"Demand_{m}"

12.10.7 Complete Transportation Example

Here’s a complete, runnable transportation problem:

from pulp import * 
 
# === DATA === 
plants = ['Chicago', 'Denver'] 
markets = ['NYC', 'LA', 'Houston'] 
 
supply = {'Chicago': 100, 'Denver': 150} 
demand = {'NYC': 80, 'LA': 90, 'Houston': 70} 
 
cost = { 
   ('Chicago', 'NYC'): 3, ('Chicago', 'LA'): 5, ('Chicago', 'Houston'): 2, 
   ('Denver', 'NYC'): 4, ('Denver', 'LA'): 2, ('Denver', 'Houston'): 3 
} 
 
# === MODEL === 
prob = LpProblem("Transportation", LpMinimize) 
 
# Variables: x[p,m] = amount shipped from plant p to market m 
routes = [(p, m) for p in plants for m in markets] 
x = LpVariable.dicts("Ship", routes, lowBound=0) 
 
# Objective: minimize total shipping cost 
prob += lpSum([cost[r] * x[r] for r in routes]), "Total_Cost" 
 
# Supply constraints 
for p in plants: 
   prob += lpSum([x[p,m] for m in markets]) <= supply[p], f"Supply_{p}" 
 
# Demand constraints 
for m in markets: 
   prob += lpSum([x[p,m] for p in plants]) >= demand[m], f"Demand_{m}" 
 
# === SOLVE === 
prob.solve() 
 
# === RESULTS === 
print(f"Status: {LpStatus[prob.status]}") 
print(f"Total Cost: ${value(prob.objective)}") 
print("\nOptimal Shipments:") 
for p in plants: 
   for m in markets: 
      if x[p,m].varValue > 0: 
         print(f" {p} -> {m}: {x[p,m].varValue}")

Info

Complete Notebooks: For fully worked examples with data files and visualizations, see:

12.11 Putting It All Together

A typical optimization workflow integrates all the tools we’ve covered:

1.
Load data with Pandas (from CSV, Excel, or databases)
2.
Explore and clean data, checking for issues
3.
Transform data into dictionaries for PuLP
4.
Build the model (variables, objective, constraints)
5.
Solve and verify status is optimal
6.
Extract results back into Pandas for analysis
7.
Visualize with Matplotlib, NetworkX, or GeoPandas
8.
Communicate findings to stakeholders
# Typical workflow skeleton 
import pandas as pd 
import matplotlib.pyplot as plt 
from pulp import * 
 
# 1. Load data 
plants_df = pd.read_csv('plants.csv') 
markets_df = pd.read_csv('markets.csv') 
costs_df = pd.read_csv('costs.csv') 
 
# 2-3. Process into dictionaries 
supply = dict(zip(plants_df['Name'], plants_df['Capacity'])) 
demand = dict(zip(markets_df['Name'], markets_df['Demand'])) 
cost = {(row['From'], row['To']): row['Cost'] 
      for _, row in costs_df.iterrows()} 
 
# 4-5. Build and solve model 
prob = LpProblem("MyProblem", LpMinimize) 
# ... define variables, objective, constraints ... 
prob.solve() 
 
# 6. Analyze results 
results_df = pd.DataFrame([...]) # Extract solution 
print(results_df.describe()) 
 
# 7. Visualize 
results_df.plot(kind='bar') 
plt.show()

The power of Python is that all these tools work seamlessly together, letting you go from raw data to optimized solution to presentation-ready visualization in a single, reproducible notebook.

Resources

Additional Resources

12.12 Exercises

Warm-ups

Exercise 12.1: A New Profit Margin

⋆  In the product mix example, a price increase raises the profit of Product 1 from $3 to $5 per unit; Product 2 still earns $2, and the resource data are unchanged. Modify the PuLP code from the example, re-solve, and report the new production plan and profit. Which constraints are binding at the new optimum? (The answer is not (20,20) anymore.)

[Section 12.10.3]

Exercise 12.2: From Math to PuLP

⋆  Translate the following linear program into PuLP code, following the six-step modeling pattern, and solve it. Report the optimal objective value and the values of all three variables.

max ⁡ 4x1 + 6x2 + 5x3  s.t.  2x1 + x2 + x3 ≤ 9 x1 + 3x2 + 2x3 ≤ 14 x1 + x2 + 2x3 ≤ 10 x1,x2,x3 ≥ 0

[Section 12.10.2, Section 12.10.3]

Exercise 12.3: New Demand Data

⋆  In the complete transportation example, market demand shifts: NYC now requires 100 units, LA requires 80, and Houston still requires 70. Supplies and costs are unchanged. Update the demand dictionary, re-solve, and report the new total cost and shipment plan. The original data gave a total cost of $610; explain in one sentence why the new cost is higher.

[Section 12.10.7]

Core problems

Exercise 12.4: A Three-Plant Transportation Model

⋆ ⋆  Build a PuLP model from scratch for the following transportation problem. Plants A, B, C can supply 60, 70, and 80 units. Markets M1, M2, M3 require 50, 60, and 70 units. Unit shipping costs are:

M1 M2 M3
A 4 6 8
B 5 4 3
C 6 7 5
Table 12.2: Unit shipping costs from plants to markets.

Use dictionaries with tuple keys for the cost data and lpSum for the objective and constraints. Report the minimum total cost, the shipment plan, and which plants have unused capacity.

[Section 12.10.7, Section 12.10.6]

Exercise 12.5: Duals and Slacks from a Solved Model

⋆ ⋆  Solve the original product mix example, then run:

for name, c in prob.constraints.items(): 
   print(name, "dual =", c.pi, "slack =", c.slack)

PuLP stores the dual value (shadow price) of each constraint in c.pi and its slack in c.slack.

1.
Report the dual value and slack of each of the three constraints.
2.
One constraint has zero dual value. What is its slack, and why does that combination make sense?
3.
If one more hour of labor became available, by how much would profit increase? Answer using only the printed output.

[Section 12.10.3, Chapter 10, Chapter 11]

Exercise 12.6: An Epsilon-Constraint Sweep

⋆ ⋆  The multi-objective chapter (Section 13.3) gives a PuLP listing that sweeps a bound on a secondary objective to trace a Pareto frontier. Use that listing as a template for the product mix example: the company wants high profit but also wants to limit machine hours 2X1 + 6X2. For each cap 𝜀 ∈{60,90,120,150,180}, maximize profit subject to the three resource constraints plus 2X1 + 6X2 ≤ 𝜀. Print a table of (𝜀,X1,X2,profit) and describe the trade-off between machine usage and profit.

[Section 13.3, Section 12.10.3]

Concepts and connections

Exercise 12.7: The Constraint That Wasn’t There

⋆ ⋆  Suppose you build the product mix model but forget to add the labor constraint (the line prob += 4*x1 + 4*x2 <= 160, "Labor" is missing).

1.
Does PuLP raise an error, print a warning, or solve without complaint? Predict, then try it.
2.
Compute the solution PuLP returns and explain why its profit is higher than $100. Is that plan actually implementable by the company?
3.
Relatedly, what happens if you type prob += 4*x1 + 4*x2 and forget the <= 160? Explain why this bug is even more dangerous than an omitted line.

[Section 12.10.3]

Exercise 12.8: Minimize or Maximize?

⋆ ⋆  A student creates the product mix problem with LpProblem("Product_Mix", LpMinimize) but keeps the objective 3*x1 + 2*x2.

1.
What solution and objective value does PuLP return? Why is the model not infeasible even though the answer is useless?
2.
Give two distinct one-line fixes: one that changes the sense and one that changes the objective. Verify in PuLP that both report the same optimal plan.
3.
Explain why checking LpStatus[prob.status] alone would not catch this bug, and suggest a sanity check that would.

[Section 12.10.2]

Challenge problems

Exercise 12.9: Profit as a Function of Labor

⋆ ⋆⋆  Write a loop that solves the product mix example for labor availability b ∈{100,110,120,… ⁡,190} (the other data unchanged), records the optimal profit for each b, and plots or tabulates profit versus b.

1.
Report the ten optimal profits. The resulting curve is piecewise linear and concave; identify the values of b where its slope changes.
2.
Print the labor constraint’s dual value c.pi at each b and explain the slopes of the curve using these shadow prices.
3.
Why does the curve become flat for large b? Answer in terms of which constraints are binding, and connect your findings to the sensitivity analysis chapter.

[Section 12.10.3, Chapter 10]

Selected Solutions

Solution

(Exercise 12.2) Following the modeling pattern:

from pulp import * 
 
prob = LpProblem("Three_Var", LpMaximize) 
x1 = LpVariable("X1", lowBound=0) 
x2 = LpVariable("X2", lowBound=0) 
x3 = LpVariable("X3", lowBound=0) 
 
prob += 4*x1 + 6*x2 + 5*x3, "Objective" 
prob += 2*x1 + x2 + x3 <= 9, "C1" 
prob += x1 + 3*x2 + 2*x3 <= 14, "C2" 
prob += x1 + x2 + 2*x3 <= 10, "C3" 
 
prob.solve() 
print(LpStatus[prob.status], value(prob.objective)) 
for v in prob.variables(): 
   print(v.name, "=", v.varValue)

Output: status Optimal, objective 35.0, with X1 = 2.0, X2 = 2.0, X3 = 3.0. All three constraints are binding ( 4 + 2 + 3 = 9, 2 + 6 + 6 = 14, and 2 + 2 + 6 = 10).

Solution

(Exercise 12.5) (1) Running the loop on the solved model prints 

Material dual = 0.2 slack = -0.0
Labor dual = 0.25 slack = -0.0
Machine dual = -0.0 slack = 20.0

(2) The machine constraint has dual value 0 and slack 20: at the optimum (20,20), machine usage is 2(20) + 6(20) = 160 < 180. A constraint with leftover capacity cannot be worth anything at the margin, so its shadow price is zero. This is complementary slackness: for each constraint, the dual value or the slack (or both) must be zero. (3) The labor dual is 0.25, so one extra hour of labor raises the optimal profit by $0.25, from $100 to $100.25 (valid while the current basis stays optimal).

Solution

(Exercise 12.9) Wrapping the model in a function of b and looping gives: 

b = 100: profit = 75.0   labor dual = 0.75
b = 110: profit = 82.5   labor dual = 0.75
b = 120: profit = 90.0   labor dual = 0.75
b = 130: profit = 92.5   labor dual = 0.25
b = 140: profit = 95.0   labor dual = 0.25
b = 150: profit = 97.5   labor dual = 0.25
b = 160: profit = 100.0  labor dual = 0.25
b = 170: profit = 102.0  labor dual = 0.0
b = 180: profit = 102.0  labor dual = 0.0
b = 190: profit = 102.0  labor dual = 0.0

(1) The slope changes at b = 120 and b = 168 (between the printed points 160 and 170). (2) The slope of the profit curve at each b equals the shadow price of the labor constraint. For b ≤ 120 only labor limits production (the plan is (b∕4,0)) and each labor hour is worth 0.75; for 120 ≤ b ≤ 168 material and labor are both binding and an extra hour is worth only 0.25; the decreasing slopes are exactly why the curve is concave. (3) At b = 168 the plan reaches (18,24), where material and machine are binding and labor is no longer scarce. Beyond that, extra labor has shadow price 0 and profit stays at $102: buying more of a non-binding resource is worthless. This is the RHS sensitivity behavior of Chapter 10, computed by brute force.

© 2026 Robert Hildebrand and contributors · Licensed CC BY-SA 4.0 · Sources and attribution · Book home