Transportation Problem (capacity allocation) in Python using Gurobi

From calculating the MODI (modified distribution) method in Operation Research (Transportation Problem) in 7 full pages during engineering to get to an optimal value, to realizing, this can be built and solved in a few minutes in excel. We have seen the evolution of mathematics with help of increased computational power.

To set a context, Transportation Problem, is a linear programming (LP) problem that identifies an optimal solution for transporting one type of product from sources (i) to destinations (j) at the minimum cost (z).

The mathematical formulation is as below:

Although this is famously known as a Transportation problem, I used this for capacity allocation when I worked as Demand & Supply Planner.

Now, the limitation with excel solver is, it can only handle 200 variables at a time, with add-ins like open solver, etc, you can increase the size of the model, but the run time will get significantly impacted.

Using python you can build the same model while there are several free solvers available in the arsenal. PuLP is the most famous free optimization modeler, there is Gurobi and even CPLEX. While a free version of the Gurobi modeler and solver license has some size limitations, the academic version comes with no size limitations.

I built the same problem using gurobi solver. Here's how I did it,

First, install the gurobi package in python:

pip install gurobipy==10.0.1        

Import the packages as you might need.

import gurobipy as gp
from gurobipy import GRB        

Defining the problem. I used limited data, to begin with, like 2 supply nodes and 4 demand nodes.

# Define data
supply = [100000, 100000]? # Supply nodes
demand = [40500, 22300, 85200, 47500]? # Demand nodes
cost = [
? ? [52, 32, 11, 69],
? ? [45, 84, 76, 15],
]? # Cost of transportation from supply i to demand j        

Create the model:

m = gp.Model()        

Define decision variables (xij), 2 represents supply nodes (i) and 4 represents demand nodes (j). lb=0 is to set the lower bound 0 considering we don't want the decision variables to be negative. Usually, in all optimization problems in the supply chain, decision variables are defined to be non-negative, however, in special cases, this can be allowed to go negative as well. e.g Inventory backorder etc.. vtype=GRB.INTEGER is to define the decision variables to be Integer (in our context), this can be floating numbers as well.

# Define decision variables
flow = m.addVars(2, 4, lb=0, vtype=GRB.INTEGER, name="flow")        

Adding supply and demand constraints. This is done simply with "For Loop" which iterates over all possible i and j as defined in the below code. Importantly, (<=) for supply constraint is added so that I do not go beyond the defined capacities or available supplies. While (>=) in demand constraint is added so that all demands are met. We can change this to (=) as well. In excel, I have observed some times, the model does anything to save the final objective value, and in order to do that, it will increase the demand in surplus if there are supplies available. gp.quicksum does the summation of all flows at destination j from supply i and ensures that LHS doesn't exceed the RHS in supply constraint while otherwise in demand constraint.

# Define supply constraints
for i in range(2):
? ? m.addConstr(gp.quicksum(flow[i, j] for j in range(4)) <= supply[i], name=f"supply_{i}")


# Define demand constraints
for j in range(4):
? ? m.addConstr(gp.quicksum(flow[i, j] for i in range(2)) >= demand[j], name=f"demand_{j}")        

Now, define the objective function. This part of the code does sumproduct of all flows from i to j, and the cost incurred doing that. sense=GRB.MINIMIZE defines that we want to minimize the cost of moving products.

# Define objective function
m.setObjective(gp.quicksum(flow[i, j] * cost[i][j] for i in range(2) for j in range(4)), sense=GRB.MINIMIZE)        

Add m.optimize to finish the model.

# Solve model
m.optimize()        

Now, printing the output needs effort, as compared to excel solver. You need to write for loop and iterate over all the flows created as output from the model, then define the printing function to print the outputs with desired notations.

# Print results
if m.status == GRB.OPTIMAL:
? ? print(f"Optimal solution found with objective value {m.objVal:.2f}")
? ? for i in range(2):
? ? ? ? for j in range(4):
? ? ? ? ? ? if flow[i, j].x > 0:
? ? ? ? ? ? ? ? print(f"Flow from supply node {i+1} to demand node {j+1}: {flow[i, j].x:.0f}")
else:
? ? print("No solution found.")        

Now, the beauty of the model is that if you want to scale it for 2000 variables you can do it, with minor modifications. If you want to supply the inputs from a CSV file, you can add them through the pandas dataframe, and tweak the codes as len(supply) and len(demand), etc. wherever needed. Change the codes for notations and it will be taken care of.

Now, an extension of the problem, where there is a processing center in between the supply and demand nodes, like Transshipment problem, needs one more additional constraint of flow conservation, and a model can be built and solved.