Data-First Optimization
Development with Gurobi andPandas
Irv Lustig, PhD, Optimization Principal, Princeton Consultants
Simon Bowly, PhD, Software Developer, Gurobi

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

1www.princeton.com

Information Technology & Management Consulting
AI & Optimization
Software Development

Strategy
Process Improvement

Princeton New Jersey

New York City

▪ Mission: Help our clients optimize their performance by integrating AI into
their core business processes▪ Track Record: 43 years in businessOver 1800 projects successfully implemented into full production acrossmultiple industriesEarned long-term strategic relationships with top transportation clients▪ Experience and Retention: Senior Staff (top 28 consultants) average 15+ years atPrinceton; Firm Leadership averages 20+ years; Client Retention > 90% annually
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

2

Dr. Irv Lustig, Optimization Principal, Princeton Consultants› Director of Numerical Optimization,
CPLEX Optimization, Inc.
› ILOG Optimization Product Manager› ILOG Optimization Evangelist› ILOG Direct Director of TechnicalServices› Vice President, ILOG Direct› Worldwide Optimization and Supply
Chain Technical Sales Manager, IBM› Manager, Optimization andMathematical Software, IBM Research

INFORMS Certified Analytics Professional
INFORMS Fellow

Core Team

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

3

The Solution Development Life Cycle› Business Problem (Question) Framing
INFORMS
CertifiedAnalytics
Professional
Job TaskAnalysis› Analytic Problem Framing› Data› Methodology (Approach) Selection› Model Building› Deployment› Life Cycle Management
See https://info.informs.org/jta for more info on the JTA
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

4

DATA FIRST MODELDEVELOPMENT
A 7-step approach that
improves the likelihood of
success in optimization
projects. This approach:• Creates useful
documentation• Shortens the overall
development time
• Improves code
maintainability• Provides a natural feedback
loop to business
requirements
7/23/2024

1. Understand the Data
2. Define and Document Datasets
3. Collect Data from the Client
4. Document Data Transformations
5. Create Model Documentation
6. Code Optimization Model

7. Validate Results and Iterate
COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

5

A Simple Example

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

6

Gurobipy Example: netflow.py› From the Gurobi Supplied Examples› A simple multicommodity network flow problem› Two products ('Pencils' and 'Pens’)– Produced in 2 cities ('Detroit' and 'Denver’)– Must be sent to warehouses in 3 cities ('Boston', 'New York', and 'Seattle’)› Satisfy supply/demand› Flows on the transportation network must respect arc capacityconstraints› The objective is to minimize the sum of the arc transportation
costs
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

7Gurobi Example: A Poor Data Representation!
# Base data
commodities = ["Pencils", "Pens"]
nodes = ["Detroit", "Denver", "Boston", "New York", "Seattle"]arcs, capacity = gp.multidict({("Detroit", "Boston"): 100,("Detroit", "New York"): 80,("Detroit", "Seattle"): 120,("Denver", "Boston"): 120,("Denver", "New York"): 120,("Denver", "Seattle"): 120,})# Cost for triplets commodity-source-destination
cost = {("Pencils", "Detroit", "Boston"): 10,("Pencils", "Detroit", "New York"): 20,("Pencils", "Detroit", "Seattle"): 60,("Pencils", "Denver", "Boston"): 40,("Pencils", "Denver", "New York"): 40,("Pencils", "Denver", "Seattle"): 30,("Pens", "Detroit", "Boston"): 20,("Pens", "Detroit", "New York"): 20,("Pens", "Detroit", "Seattle"): 80,("Pens", "Denver", "Boston"): 60,("Pens", "Denver", "New York"): 70,("Pens", "Denver", "Seattle"): 30,}
7/23/2024# Supply (> 0) and demand (< 0) for pairs ofcommodity-city
inflow = {("Pencils", "Detroit"): 50,("Pencils", "Denver"): 60,("Pencils", "Boston"): -50,("Pencils", "New York"): -50,("Pencils", "Seattle"): -10,("Pens", "Detroit"): 60,("Pens", "Denver"): 40,("Pens", "Boston"): -40,("Pens", "New York"): -30,("Pens", "Seattle"): -30,}

What’s Poor About This?
• Data inside code
• Dependent on Gurobi multidict• Commodities and Nodes repeated• These can be inferred!
• Mixing Supply and Demand Data
• If there were transshipment nodes,
inflow could have many zero values
COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

8

A Better Data Representation: Tables

1. Understand the Data

3. Collect Data from the Client

import pandas as pd
raw_data = pd.read_excel("mcfdata.v1.xlsx", sheet_name=None)
data = {k: df.set_index(df.columns[:-1].to_list())[df.columns[-1]]for k, df in raw_data.items()}
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

9

Pandas Representation: A Dictionary of Series

7/23/2024

1. Understand the Data

3. Collect Data from the Client

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

10Document the Data Model

2. Define and Document Datasets

Sets
Notation
𝒩Meaning
Set of network nodes𝒫

Set of products (commodities)

𝒜𝒫𝑎
𝒜𝑝

Set of arcs (𝑛𝑓 , 𝑛𝑡 ), 𝑛𝑓 , 𝑛𝑡 ∈ 𝒜
Set of products 𝑝 ∈ 𝒫 that can be carried on arc 𝑎 ∈ 𝒜
Set of arcs 𝑎 ∈ 𝒜 that can carry product 𝑝 ∈ 𝒫

Table Locations
cost: Columns From, To
capacity: Columns From, To
supply: Column Nodedemand: Column Nodecost: Column Productsupply: Column Productdemand: Column Productcost: Columns From, To
cost: Columns Product, From, To
cost: Columns Product, From, To

Numerical Input Values
The input data is converted to pandas Series, so the name of each Series is also the name of the value.

Notation Meaning𝜅𝑎Capacity of arc 𝑎 ∈ 𝒜
𝜋𝑎𝑝
Cost of carrying product 𝑝 on arc 𝑎 ∈ 𝒜, 𝑝 ∈ 𝒫𝑎 ,𝜎𝑝𝑛
𝛿𝑝𝑛

Supply of product 𝑝 ∈ 𝒫 at node 𝑛 ∈ 𝒩. Defaults
to 0
Demand of product 𝑝 ∈ 𝒫 at node 𝑛 ∈ 𝒩.
Defaults to 07/23/2024

Table Name/Value Columncapacity
cost

Index Columns
From, To
Product, From, To

supplyNodedemand

Node

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

11Computing Sets from Data

4. Document Data Transformations› In a Data First approach, sets are inferred from data
– The set 𝒫 of products can appear in any of the tables: supply, demand and cost .
– The set 𝒩 of nodes can appear in any of the tables: capacity, supply, demand and costcommodities = set(pd.concat([data[dfname].index.to_frame()["Product"]
for dfname in ["supply", "demand", "cost"]]).unique())
nodes = set(pd.concat([data[dfname].index.to_frame()[fromto].rename("Node")for dfname in ["capacity", "cost"]
for fromto in ["From", "To"]]+ [data[dfname].index.to_frame()["Node"] for dfname in ["supply", "demand"]]).unique())
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

12

Computing the Net Flow Per Node

4. Document Data Transformations› The net flow 𝜇𝑝𝑛 for each product 𝑝 ∈ 𝒫 and node 𝑛 ∈ 𝒩 is the sum of the supply lessthe demand. For transshipment nodes, this value is 0. This is called inflow in thecode.
inflow = pd.concat([data["supply"].rename("net"),data["demand"].rename("net") * -1,pd.Series(
0,index=pd.MultiIndex.from_product([commodities, nodes], names=["Product", "Node"]),name="net",),]).groupby(["Product", "Node"]).sum()› Note that the above allows supply and demand to co-exist at any node, andfor any node to act as transshipment nodes› This is not true in the original Gurobi example.
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

13

Time to Model

5. Create Model Documentation› Decision Variables–

The model will have one set of decision variables: 𝑋𝑝𝑎 for 𝑝 ∈ 𝒫, 𝑎 ∈ 𝒜𝑝 represents the amount shipped of product 𝑝 on arc 𝑎. We
will call this variable flow in the code.

–

The cost of shipment is 𝜋𝑎𝑝 .

–This defines the objective function:minimize෍ ෍ 𝜋𝑎𝑝 𝑋𝑝𝑎𝑎∈𝒜 𝑝∈𝒫𝑎

› Constraints–

Flow on each arc is capacitated෍ 𝑋𝑝𝑎 ≤ 𝜅𝑎∀𝑎 ∈ 𝒜𝑝∈𝒫𝑎–Conservation of Flow–

For each node and each product, the flow out of the node, less the flow into the node is equal to the net flow.

෍(𝑛,𝑛𝑡 )∈𝐴𝑝
7/23/2024

𝑋𝑝(𝑛,𝑛𝑡) − ෍ 𝑋𝑝(𝑛𝑓,𝑛) = 𝜇𝑝𝑛∀𝑝 ∈ 𝒫, 𝑛 ∈ 𝒩(𝑛𝑓 ,𝑛)COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

14

Using gurobipy-pandas: Decision Variables

6. Code Optimization Model

import gurobipy as grb
import gurobipy_pandas as gppd
m = grb.Model("netflow")
flow = gppd.add_vars(m, data["cost"], obj=data["cost"], name="flow")
Decision Variables ShareIndex with pandas Series

ProductPencils

FromDetroitToBoston
New YorkSeattle
Denver
Boston
New York
SeattlePens
Detroit Boston
New YorkSeattle
Denver
Boston
New York
SeattleName: flow, dtype: object
7/23/2024

Series has the values forthe objective<gurobi.Var flow[Pencils,Detroit,Boston]><gurobi.Var flow[Pencils,Detroit,New_York]><gurobi.Var flow[Pencils,Detroit,Seattle]><gurobi.Var flow[Pencils,Denver,Boston]><gurobi.Var flow[Pencils,Denver,New_York]><gurobi.Var flow[Pencils,Denver,Seattle]><gurobi.Var flow[Pens,Detroit,Boston]><gurobi.Var flow[Pens,Detroit,New_York]><gurobi.Var flow[Pens,Detroit,Seattle]><gurobi.Var flow[Pens,Denver,Boston]><gurobi.Var flow[Pens,Denver,New_York]><gurobi.Var flow[Pens,Denver,Seattle]>
COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

15

Using gurobipy-pandas: Capacity Constraints

6. Code Optimization Model෍ 𝑋𝑝𝑎 ≤ 𝜅𝑎∀𝑎 ∈ 𝒜𝑝∈𝒫𝑎

pd.concat([flow.groupby(["From", "To"]).agg(grb.quicksum), data["capacity"]], axis=1).gppd.add_constrs(m, "flow <= capacity", name="cap")7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

16

Using gurobipy-pandas: Flow Constraints

6. Code Optimization Model෍(𝑛,𝑛𝑡 )∈𝐴𝑝

𝑋𝑝(𝑛,𝑛𝑡 ) − ෍ 𝑋𝑝(𝑛𝑓 ,𝑛) = 𝜇𝑝𝑛(𝑛𝑓 ,𝑛)pd.concat([flow.rename_axis(index={"From": "Node"}).groupby(["Product", "Node"]).agg(grb.quicksum).rename("flowout"),

flow.rename_axis(index={"To": "Node"}).groupby(["Product", "Node"]).agg(grb.quicksum).rename("flowin"),inflow,∀𝑝 ∈ 𝒫, 𝑛 ∈ 𝒩෍

𝑋𝑝(𝑛,𝑛𝑡 )(𝑛,𝑛𝑡 )∈𝐴𝑝

෍ 𝑋𝑝(𝑛𝑓 ,𝑛)(𝑛𝑓 ,𝑛)𝜇𝑝𝑛],
axis=1,).fillna(0).gppd.add_constrs(m, "flowout - flowin == net", name="node")

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

17

Flow Constraint in pandas

6. Code Optimization Model› Note the 0 values are created via fillna(0) because there isno flow into supply nodes and flow out of demand nodes, so thepd.concat operation would create NA values, which cannot bepassed to Gurobi
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

18

Optimize and Get Results

7. Validate Results and Iteratem.optimize()soln = flow.gppd.X
soln.to_frame().query("flow > 0").sort_index()› Note the use of the gppd accessor to get the attribute X of the decisionvariables flow use– All Gurobi attributes can be accessed in this way!

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

19

Energy Systems Modeling: Unit Commitment Problem› Dispatch generators to meet expected demand› Keep a target level of available capacity online› Minimize overall generation cost
› Many more complex constraints in real life!
gurobipy-pandas.readthedocs.io

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

20

Pandas Representation: DataFrames› Generators have:– Fixed costs and marginal costs– Startup costs– Lower and upper generation limits when online› Each time period has:
– Demand to be met– A minimum online capacity› Model captures state of each class in each time period
7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

21

Decision Variables: gurobipy-pandas
generators = (# Decision variables for each class at each time period
pd.DataFrame(index=pd.MultiIndex.from_product([generator_data.index,time_period_data.index,])).join(generator_data)# Three classes of variables capture output and on/off state.gppd.add_vars(model, name="output").gppd.add_vars(model,vtype=GRB.INTEGER,
ub="num_available",
name="num_active",).gppd.add_vars(model, vtype=GRB.INTEGER, name="num_startup”))

7/23/2024

Combined timeexpanded index(with generatordata)

Each accessor callcreates one
variable per row

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

22

Time-Expanded Format
Time-expanded
pandas MultiIndex

Decision variables share
this common index

7/23/2024

Generator data is
expanded (makes modelbuilding easier)

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

23

Constraints: Generator Limits
Time-expanded
dataframe# Generators must operate within their limits
generators.gppd.add_constrs(model,"output >= min_output * num_active",
name="lower_limit",).gppd.add_constrs(
model,"output <= max_output * num_active",
name="upper_limit",)

7/23/2024

One constraint
per timeexpanded row

New constraint defines
a relationship between
columns

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

24

Constraints: Satisfying Demand# Predicted demand must be exactly satisfied
demand_constraint = gppd.add_constrs(
model,generators.groupby("time_period")["output"].sum(),
GRB.EQUAL,
time_period_data["expected_demand"],)

These two
expressions havecompatible indexes

groupby + sum createsa constraint acrossgenerators for one
time period

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

25

Constraints: Reserve Capacity# Ensure that the active generators during each time
# period are capable of meeting the reserve demand.
active_capacity = ((generators["max_output"] * generators["num_active"]).groupby("time_period").sum())
active_capacity_constraint = gppd.add_constrs(
model,
active_capacity,
GRB.GREATER_EQUAL,
time_period_data["minimum_active_capacity"],)Keep this forresults step

Expression
measuring total
active capacity in agiven time period

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

26

Constraints: Startup Costs
def count_startups(group):
return gppd.add_constrs(
model,group["num_startup"].iloc[1:],GRB.GREATER_EQUAL,group["num_active"].diff().dropna(),
name="startup",).to_frame()diff() captures thechange from one time
period to the next

startup = generators.groupby("generator_class").apply(count_startups)

groupby creates a
set of constraints
per generator class

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

27

Objective: Total Cost
# Minimize total cost objectivemodel.setObjective((# Fixed hourly costs for started generators
generators["cost_per_hour"] * generators["num_active"]# Marginal hourly cost of generation above the minimum+ generators["marginal_cost"] * (generators["output"] –generators["num_active"] * generators["min_output"])# Startup costs for newly active generators in each period+ generators["startup_cost"] * generators["num_startup"]).sum(),
sense=GRB.MINIMIZE
Aggregate)
objective

7/23/2024

Creates a cost
expression perperiod and class

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

28

Optimize and Extract Resultsmodel.optimize()
results = pd.DataFrame({"Demand": time_period_data["expected_demand"],"Min. Active Capacity": time_period_data["minimum_active_capacity"],"Active Capacity": active_capacity.gppd.get_value(),"Excess": (-active_capacity_constraint.gppd.Slack),})

Value ofexpressions inthe solutionresults.plot.line()

Slack in constraints(Gurobi attribute query)

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)

29

CONCLUSION

1. Understand the Data

• Data First Model
Development Focuseson the Data

2. Define and Document Datasets• A proven practice for delivering highquality optimization applications• gurobipy-pandas• Takes advantage of pandas for
managing data
• Allow Gurobi objects to be placed intopandas DataFrames and Series• Provides capability to write code that
executes quickly when creatingmodels• Requires good knowledge of pandas
to understand how to map
mathematical formulations into
gurobipy-pandas
7/23/2024

3. Collect Data from the Client
4. Document Data Transformations
5. Create Model Documentation
6. Code Optimization Model

7. Validate Results and Iterate
COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)30

Thanks!
› More info:
o https://gurobipy-pandas.readthedocs.io/en/stable/o https://www.princetonoptimization.com/blog/blog/rapid-optimizationmodel-development-python-and-pandas-7-steps

7/23/2024

COPYRIGHT © PRINCETON CONSULTANTS AND GUROBI 2024 (V1.0)31