Home About
Courses Concepts Tools
GitHub

Min Cost Flow ILP Formulation

Min Cost Flow ILP (GurobiPy)

Author: Ben Rosenberg

Imports

We begin by importing some relevant libraries. GurobiPy we import to formulate and solve the ILP; time, we import to time the entire process of supplying the constraints solving the ILP.

import gurobipy as gp
from gurobipy import GRB
import time

Input data

Next, we define our input data. Recall that in the Min Cost Flow problem we have as an input a graph G = (N, E), where N is the set of nodes, each of which has supply s(i) for all i\in N, and E is the set of edges, each of which has a cost c(i,j) and capacity u(i,j) for all (i,j)\in E

The input data below corresponds to a small example graph, seen below, in which the numbers inside the nodes are their supply/demand values, and the first element of the tuple corresponding to each edge denotes the cost and the second denotes the capacity:

nodes = range(4)

edges = [(0,1), (0,3), (2,1), (2,3)]

# cost[i,j] is the cost of sending one unit of flow along edge (i,j)
cost = {
    (0,1) : 0,
    (0,3) : 4,
    (2,1) : 2,
    (2,3) : 8
} 

# capacity[i,j] is the capacity of edge (i,j)
capacity = {
    (0,1) : 5,
    (0,3) : 2,
    (2,1) : 3,
    (2,3) : 2
}

# supply[i] is the supply (if positive) or demand (if negative) of node i
supply = [5, -6, 3, -2]

Model Definition

Now we define our model. The Min Cost Flow problem has the following constraints, using the variable f(i,j) to denote the amount of flow on edge (i,j):

  • Flow in is equal to flow out, adjusting for the offset of supply \sum_{(i,j)\in E} f(i,j) - \sum_{(j,i)\in E} f(j,i) = s(i) \quad \forall i\in N
  • Flow on an edge can be no more than the capacity of that edge f(i,j) \leq u(i,j) \quad \forall (i,j)\in E
  • Flow on an edge can be no less than 0 units f(i,j) \geq 0 \quad \forall (i,j) \in E

And the objective should be intuitive, as it's simply minimizing the cost of all flows on the edges:

\min \sum_{(i,j)\in E} c(i,j) \cdot f(i,j)

start_time = time.time()

model = gp.Model()

# set output level to max
model.Params.TuneOutput = 3

# add variable f
f = model.addVars(nodes, nodes, vtype=GRB.CONTINUOUS, name='f')

# add constraint representing supply/demand
for i in nodes:
    model.addConstr(gp.quicksum(f[i,j] for (x,j) in edges if x == i) -
                    gp.quicksum(f[j,i] for (j,x) in edges if x == i)
                    == supply[i])

# add constraint on edge flows w.r.t. capacities
for (i,j) in edges:
    model.addConstr(f[i,j] <= capacity[i,j])

# add constraint on edge flows w.r.t. 0
for (i,j) in edges:
    model.addConstr(f[i,j] >= 0)
    
# set objective
model.setObjective(gp.quicksum(cost[i,j] * f[i,j] for (i,j) in edges), GRB.MINIMIZE)

model.optimize()

end_time = time.time()

diff = time.gmtime(end_time - start_time)
print('\n[Total time used: {} minutes, {} seconds]'.format(diff.tm_min, diff.tm_sec))

try:
    print(f'\nObjective value found: {model.objVal}')
except AttributeError as e:
    print(f'\nCould not find an objective value. \nTraceback:\n\t{e}')

Set parameter Username
Academic license - for non-commercial use only - expires 2023-10-13
Set parameter TuneOutput to value 3
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 12 rows, 16 columns and 16 nonzeros
Model fingerprint: 0xaf2c58bd
Coefficient statistics:
    Matrix range     [1e+00, 1e+00]
    Objective range  [2e+00, 8e+00]
    Bounds range     [0e+00, 0e+00]
    RHS range        [2e+00, 6e+00]
Presolve removed 12 rows and 16 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
        0    1.4000000e+01   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.01 seconds (0.00 work units)
Optimal objective  1.400000000e+01

[Total time used: 0 minutes, 0 seconds]

Objective value found: 14.0

So we have our optimal objective. The optimal solution associated with that, in terms of flows on edges, is given below:

# optimal solution
for (i,j) in edges:
    print('f[{},{}] = {}'.format(
        i, j, f[i,j].x
    ))

f[0,1] = 3.0
f[0,3] = 2.0
f[2,1] = 3.0
f[2,3] = 0.0

Min Cost Flow ILP (OR-Tools)

Author: Ben Rosenberg

Imports

We begin by importing some relevant libraries. We import ortools to formulate and solve the ILP, and we import time to time the entire process of supplying the constraints solving the ILP.

from ortools.linear_solver import pywraplp as OR
import time

Input data

Next, we define our input data. Recall that in the Min Cost Flow problem we have as an input a graph G = (N, E), where N is the set of nodes, each of which has supply s(i) for all i\in N, and E is the set of edges, each of which has a cost c(i,j) and capacity u(i,j) for all (i,j)\in E

The input data below corresponds to a small example graph, seen below, in which the numbers inside the nodes are their supply/demand values, and the first element of the tuple corresponding to each edge denotes the cost and the second denotes the capacity:

nodes = range(4)

edges = [(0,1), (0,3), (2,1), (2,3)]

# cost[i,j] is the cost of sending one unit of flow along edge (i,j)
cost = {
    (0,1) : 0,
    (0,3) : 4,
    (2,1) : 2,
    (2,3) : 8
} 

# capacity[i,j] is the capacity of edge (i,j)
capacity = {
    (0,1) : 5,
    (0,3) : 2,
    (2,1) : 3,
    (2,3) : 2
}

# supply[i] is the supply (if positive) or demand (if negative) of node i
supply = [5, -6, 3, -2]

Model Definition

Now we define our model. The Min Cost Flow problem has the following constraints, using the variable f(i,j) to denote the amount of flow on edge (i,j):

  • Flow in is equal to flow out, adjusting for the offset of supply \sum_{(i,j)\in E} f(i,j) - \sum_{(j,i)\in E} f(j,i) = s(i) \quad \forall i\in N
  • Flow on an edge can be no more than the capacity of that edge f(i,j) \leq u(i,j) \quad \forall (i,j)\in E
  • Flow on an edge can be no less than 0 units f(i,j) \geq 0 \quad \forall (i,j) \in E

And the objective should be intuitive, as it's simply minimizing the cost of all flows on the edges:

\min \sum_{(i,j)\in E} c(i,j) \cdot f(i,j)

start_time = time.time()

m = OR.Solver('MCF', OR.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# add variable f
f = {}
for (i,j) in edges:
    f[i,j] = m.NumVar(0, m.infinity(), 'f[{},{}]'.format(
        i,j
    ))

# add constraint representing supply/demand
for i in nodes:
    m.Add(sum(f[x,j] for (x,j) in edges if x == i) - 
            sum(f[j,x] for (j,x) in edges if x == i)
            == supply[i])

# add constraint on edge flows w.r.t. capacities
for (i,j) in edges:
    m.Add(f[i,j] <= capacity[i,j])

# add constraint on edge flows w.r.t. 0
for (i,j) in edges:
    m.Add(f[i,j] >= 0)
    
# set objective
m.Minimize(sum(cost[i,j] * f[i,j] for (i,j) in edges))

m.Solve()

end_time = time.time()

diff = time.gmtime(end_time - start_time)
print('\n[Total time used: {} minutes, {} seconds]'.format(diff.tm_min, diff.tm_sec))

print('Objective:', m.Objective().Value())


[Total time used: 0 minutes, 0 seconds]
Objective: 14.0

So we have our optimal objective. The optimal solution associated with that, in terms of flows on edges, is given below:

# optimal solution
for (i,j) in edges:
    print('f[{},{}] = {}'.format(
        i, j, f[i,j].solution_value()
    ))

f[0,1] = 3.0
f[0,3] = 2.0
f[2,1] = 3.0
f[2,3] = 0.0