Max Flow ILP Formulation
Max 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 Max Flow problem we have as an input a graph G = (N, E), where N is the set of nodes, and E is the set of edges, each of which has a 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 arbitrary indices, and the number corresponding to each edge denotes the capacity of that edge:
nodes = range(6)
edges = [(0,1), (0,3), (1,2), (1,3), (2,3), (2,5), (3,4), (4,5)]
# capacity[i,j] is the capacity of edge (i,j)
capacity = {
(0,1) : 5,
(0,3) : 3,
(1,2) : 3,
(1,3) : 1,
(2,3) : 2,
(2,5) : 6,
(3,4) : 4,
(4,5) : 3
}
Model Definition
Now we define our model. The Max 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 \sum_{(i,j)\in E} f(i,j) - \sum_{(j,i)\in E} f(j,i) = 0 \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 maximizing the total flow through the graph (which is equivalent to maximizing the flow into the sink node, 5):
\max \sum_{(i,5)\in E} f(i,5)
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 flow in/out
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)
== 0)
# 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(f[i,j] for (i,j) in edges), GRB.MAXIMIZE)
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 22 rows, 36 columns and 32 nonzeros
Model fingerprint: 0x54b6a1d8
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e+00, 1e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 6e+00]
Presolve removed 22 rows and 36 columns
Presolve time: 0.01s
Presolve: All rows and columns removed
Iteration Objective Primal Inf. Dual Inf. Time
0 -0.0000000e+00 0.000000e+00 0.000000e+00 0s
Solved in 0 iterations and 0.01 seconds (0.00 work units)
Optimal objective -0.000000000e+00
[Total time used: 0 minutes, 0 seconds]
Objective value found: -0.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] = 0.0 f[0,3] = 0.0 f[1,2] = 0.0 f[1,3] = 0.0 f[2,3] = 0.0 f[2,5] = 0.0 f[3,4] = 0.0 f[4,5] = 0.0
Max 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 Max Flow problem we have as an input a graph G = (N, E), where N is the set of nodes, and E is the set of edges, each of which has a 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 arbitrary indices, and the number corresponding to each edge denotes the capacity of that edge:
nodes = range(6)
edges = [(0,1), (0,3), (1,2), (1,3), (2,3), (2,5), (3,4), (4,5)]
# capacity[i,j] is the capacity of edge (i,j)
capacity = {
(0,1) : 5,
(0,3) : 3,
(1,2) : 3,
(1,3) : 1,
(2,3) : 2,
(2,5) : 6,
(3,4) : 4,
(4,5) : 3
}
Model Definition
Now we define our model. The Max 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 \sum_{(i,j)\in E} f(i,j) - \sum_{(j,i)\in E} f(j,i) = 0 \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 maximizing the total flow through the graph (which is equivalent to maximizing the flow into the sink node, 5):
\max \sum_{(i,5)\in E} f(i,5)
start_time = time.time()
m = OR.Solver('Max Flow', 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 flow in/out for non-source/sink nodes
for i in nodes:
if i not in (0,5):
m.Add(sum(f[i,j] for (x,j) in edges if x == i) -
sum(f[j,i] for (j,x) in edges if x == i)
== 0)
# 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.Maximize(sum(f[i,j] for (i,j) in edges if j == 5))
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: 6.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] = 3.0 f[1,2] = 3.0 f[1,3] = 0.0 f[2,3] = 0.0 f[2,5] = 3.0 f[3,4] = 3.0 f[4,5] = 3.0