The goal of the maximum flow problem is to find the maximum flow possible in a network from some given source node to a given sink node. Applications of the max flow problem include finding the maximum flow of orders through a job shop, the maximum flow of water through a storm sewer system, and the maximum flow of product through a product distribution system, among others. Schrijver (2002) note that the maximum flow problem was first formulated in 1954 by T. E. Harris and F. S. Ross as a simplified model of Soviet railway traffic flow.
A network is a directed graph, and the arc capacities, or upper bounds, are the only relevant parameters. A network graph does not have to be symmetric: if an arc (v,w) is in the graph, the reverse arc (w,v) does not have to be in the graph. Further, parallel arcs are not allowed, and self-loops are not allowed. The source and the sink are distinct nodes in the network, but the sink may be unreachable from the source.
The max flow problem can be formulated mathematically as a linear programming problem using the following model.
$N$ = nodes in the network
$A$ = network arcs
$s$ = source node
$t$ = sink node
$c_{ij}$ = arc flow capacity, $\forall (i,j) \in A$
$f_{i,j}$ = arc flow, $\forall (i,j) \in A$
Maximize the flow into the sink nodes
$\max \sum_{\{i \mid (i,t) \in A\}} c_{i,t} f_{i,t}$
Enforce an upper limit on the flow across each arc
$f_{i,j} \leq c_{i,j}$, $\forall (i,j) \in A$
Enforce flow through each node
$\sum_{\{i \mid (i,j) \in A\}} f_{i,j} = \sum_{\{i \mid (j,i) \in A\}} f_{j,i}$, $\forall j \in N - \{s,t\}$
Flow lower bound
$f_{i,j} \geq 0$, $\forall (i,j) \in A$
We begin by importing the Pyomo package and creating a model object:
from pyomo.environ import *
model = AbstractModel()
The sets $N$ and $A$ are declared abstractly using the Set
component:
# Nodes in the network
model.N = Set()
# Network arcs
model.A = Set(within=model.N*model.N)
Similarly, the model parameters are defined abstractly using the Param
component:
# Source node
model.s = Param(within=model.N)
# Sink node
model.t = Param(within=model.N)
# Flow capacity limits
model.c = Param(model.A)
The within
option is used in these parameter declarations to define expected properties of the parameters. This information is used to perform error checks on the data that is used to initialize the parameter components.
The Var
component is used to define the decision variables:
# The flow over each arc
model.f = Var(model.A, within=NonNegativeReals)
The within
option is used to restrict the domain of the decision variables to the non-negative reals. This eliminates the need for explicit bound constraints for variables.
The Objective
component is used to define the cost objective. This component uses a rule function to construct the objective expression:
# Maximize the flow into the sink nodes
def total_rule(model):
return sum(model.f[i,j] for (i, j) in model.A if j == value(model.t))
model.total = Objective(rule=total_rule, sense=maximize)
Similarly, rule functions are used to define constraint expressions in the Constraint
component:
# Enforce an upper limit on the flow across each arc
def limit_rule(model, i, j):
return model.f[i,j] <= model.c[i, j]
model.limit = Constraint(model.A, rule=limit_rule)
# Enforce flow through each node
def flow_rule(model, k):
if k == value(model.s) or k == value(model.t):
return Constraint.Skip
inFlow = sum(model.f[i,j] for (i,j) in model.A if j == k)
outFlow = sum(model.f[i,j] for (i,j) in model.A if i == k)
return inFlow == outFlow
model.flow = Constraint(model.N, rule=flow_rule)
Putting these declarations all together gives the following model:
!cat maxflow.py
from pyomo.environ import * model = AbstractModel() # Nodes in the network model.N = Set() # Network arcs model.A = Set(within=model.N*model.N) # Source node model.s = Param(within=model.N) # Sink node model.t = Param(within=model.N) # Flow capacity limits model.c = Param(model.A) # The flow over each arc model.f = Var(model.A, within=NonNegativeReals) # Maximize the flow into the sink nodes def total_rule(model): return sum(model.f[i,j] for (i, j) in model.A if j == value(model.t)) model.total = Objective(rule=total_rule, sense=maximize) # Enforce an upper limit on the flow across each arc def limit_rule(model, i, j): return model.f[i,j] <= model.c[i, j] model.limit = Constraint(model.A, rule=limit_rule) # Enforce flow through each node def flow_rule(model, k): if k == value(model.s) or k == value(model.t): return Constraint.Skip inFlow = sum(model.f[i,j] for (i,j) in model.A if j == k) outFlow = sum(model.f[i,j] for (i,j) in model.A if i == k) return inFlow == outFlow model.flow = Constraint(model.N, rule=flow_rule)
Since this is an abstract Pyomo model, the set and parameter values need to be provided to initialize the model. The following data command file provides a synthetic data set:
!cat maxflow.dat
set N := Zoo A B C D E Home; set A := (Zoo,A) (Zoo,B) (A,C) (A,D) (B,A) (B,C) (C,D) (C,E) (D,E) (D,Home) (E,Home); param s := Zoo; param t := Home; param: c := Zoo A 11 Zoo B 8 A C 5 A D 8 B A 4 B C 3 C D 2 C E 4 D E 5 D Home 8 E Home 6;
Set data is defined with the set
command, and parameter data is defined with the param
command.
This data set considers the problem of maximizing flow in a zoo. A panda is about to give birth at the zoo! Officials anticipate that attendance will skyrocket to see the new, adorable baby panda. There's one particular residential area called "Home" that is full of panda loving families and there's a fear that the increased number of people visiting the zoo will overload the public transportation system. It will be especially bad in the evening since the zoo closes about the same time as rush hour, so everyone will be trying to find space on the already crowded buses and subways. As a city planner, you were given a map of routes from the zoo to Home, along with the estimated number of families that could go on each route. Additionally, it was estimated that 16 families from Home will visit each day, and it's your task to figure out if this will overload the public transportation system, and, if it does, how could the system be improved?
Pyomo includes a pyomo
command that automates the construction and optimization of models. The GLPK solver can be used in this simple example:
!pyomo solve --solver=glpk maxflow.py maxflow.dat
[ 0.00] Setting up Pyomo environment [ 0.00] Applying Pyomo preprocessing actions [ 0.00] Creating model [ 0.02] Applying solver [ 0.06] Processing results Number of solutions: 1 Solution Information Gap: 0.0 Status: feasible Function Value: 14.0 Solver results file: results.json [ 0.06] Applying Pyomo postprocessing actions [ 0.06] Pyomo Finished
By default, the optimization results are stored in the file results.yml
:
!cat results.yml
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 14.0 Upper bound: 14.0 Number of objectives: 1 Number of constraints: 17 Number of variables: 12 Number of nonzeros: 30 Sense: maximize # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: ok Termination condition: optimal Statistics: Branch and bound: Number of bounded subproblems: 0 Number of created subproblems: 0 Error rc: 0 Time: 0.00943398475647 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 1 number of solutions displayed: 1 - Gap: 0.0 Status: feasible Message: None Objective: total: Value: 14 Variable: f[C,E]: Value: 4 f[C,D]: Value: 2 f[Zoo,A]: Value: 7 f[A,C]: Value: 3 f[E,Home]: Value: 6 f[B,A]: Value: 4 f[D,Home]: Value: 8 f[D,E]: Value: 2 f[B,C]: Value: 3 f[Zoo,B]: Value: 7 f[A,D]: Value: 8 Constraint: No values
This output tells us how many people should travel along each route for the optimal solution. More importantly, though, is the line which says our objective value is 14. This means that at most 14 families can arrive at Home. However, we were told 16 families from Home were expected to visit the zoo each day. Therefore, unless something is done, the public transportation network in place will be overloaded.