#!/usr/bin/env python
# coding: utf-8

# # Max Flow
# 
# ## Summary
# 
# 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.
#  
# 
# ## Problem Statement
# 
# The max flow problem can be formulated mathematically as a linear programming problem using the following model.  
# 
# ### Sets
# 
#  $N$ = nodes in the network  
#  $A$ = network arcs
# 
# ### Parameters
# 
#  $s$ = source node  
#  $t$ = sink node  
#  $c_{ij}$ = arc flow capacity, $\forall (i,j) \in A$
#  
# ### Variables
#  $f_{i,j}$ = arc flow, $\forall (i,j) \in A$
# 
# ### Objective
# 
# Maximize the flow into the sink nodes  
#  $\max \sum_{\{i \mid (i,t) \in A\}} c_{i,t} f_{i,t}$
# 
# ### Constraints
# 
# 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$
# 
# ## Pyomo Formulation
# 
# We begin by importing the Pyomo package and creating a model object:

# In[1]:


from pyomo.environ import *

model = AbstractModel()


# The sets $N$ and $A$ are declared abstractly using the `Set` component:

# In[2]:


# 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:

# In[3]:


# 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:

# In[4]:


# 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:

# In[5]:


# 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:

# In[6]:


# 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:

# In[7]:


get_ipython().system('cat maxflow.py')


# ## Model Data
# 
# 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:

# In[8]:


get_ipython().system('cat maxflow.dat')


# 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?
# 
# ## Solution
# 
# Pyomo includes a `pyomo` command that automates the construction and optimization of models.  The GLPK solver can be used in this simple example:

# In[9]:


get_ipython().system('pyomo solve --solver=glpk maxflow.py maxflow.dat')


# By default, the optimization results are stored in the file `results.yml`:

# In[10]:


get_ipython().system('cat results.yml')


# 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.

# ## References
# 
# * A. Schrijver, (2002). "On the history of the transportation and maximum flow problems". Mathematical Programming 91 (3): 437–445.