```python
import gurobipy as gp
from gurobipy import GRB
```

# Max-Flow Min-Cut



The maximum flow minimum cut problem holds a special place in the history of optimization theory. We will first model the problems as Linear Programs and use the results to discuss some somewhat surprising results.

The problems consider a directed graph consisting of a set of nodes and a set of labeled arcs. The arc labels are non-negative values representing a notion of capacity for the arc. In the node set, there exists a source node $s$ and a terminal node $t$. The amount of flow into one of the intermediary nodes must equal the amount of flow out of the node, i.e., flow is conserved. 

The maximum flow question asks: what is the maximum flow that can be transferred from the source to the sink. The minimum cut question asks: which is the subset of arcs, that once removed would disconnect the source node from the terminal node, which has the minimum sum of capacities. For example, removing arcs $(C, t)$ and $(D, t)$ from the network below would mean there is no longer a path from $s$ to $t$ and the sum of the capacities of these arcs is $140 + 90 = 250$. It is reasonably straight forward to find a better cut, i.e., a subset of nodes with sum of capacities less than 250.

A complete model is provided for the maximum flow problem, whereas the minimum cut problem is left as a challenge to the reader.




### Notation

Let's represent the set of nodes and arcs as below

|  index   |  Set   |  Description | 
|:---------|:--------|:--------------|
| $i$ | $V$ | Set of nodes  |
| $(i,j)$ |  $A$   | Set of arcs  |

Additionally we can define the capacity of the arcs as follows

|  Parameter |  Description | 
|:---------|:--------|
| $c_{i,j}$ | Capacity of arc $(i,j) \in A$  |


```python
# Programmatically we can define the problem data as 
nodes = ['s', 'A', 'B', 'C', 'D', 't']
capacity = {
    ('s', 'A'):   100,
    ('s', 'B'):  150,
    ('A', 'B'):  120,
    ('A',  'C'):   90,
    ('B',  'D'): 110,
    ('C',  'D'):  120,
    ('C', 't'): 140,
    ('D', 't'): 90,
    }
arcs = capacity.keys()
```

## Maximum Flow

First, let's consider the problem of calculating the maximum flow that can pass through the network.

### Variables

The variables we will use are as follows

|  Variable | Type | Description | 
|:---------|:--------| :----- |
| $f_{i,j}$ | Continuous | Flow along arc $(i,j) \in A$  |

### Model

A model of the problem can then be defined as follows:

$$
\begin{align}
\text{maximise} \ & \sum_{j \in V: (s,j) \in A} f_{s, j} && & \quad (1a) \label{model-obj}\\
s.t. \ & f_{i, j} \leq c_{i,j} \quad && \forall (i, j) \in A & \quad (1b) \label{m1-c1}\\
& \sum_{i \in V: (i,j) \in A} f_{i, j} - \sum_{k \in V: (j, k) \in A} f_{j,k} = 0 \quad  && \forall j \in V \setminus \{s, t\} & \quad (1c) \label{m2-c2}
\end{align}
$$

The objective (1a) is to maximise the sum of flow leaving the source node $s$. Constraints (1b) ensure that the flow in each arc does not exceed the capacity of that arc. Constraints (1c) are continuity constraints, which ensure that the flow into each of the nodes, excluding the source and sink, is equal to the flow out of that node.


```python
# A function that takes a set of nodes, arcs, and capacities, creates a model and optimises can be defined as follows:
def max_flow(nodes, arcs, capacity):

    # Create optimization model
    m = gp.Model('flow')

    # Create variables
    flow = m.addVars(arcs, obj=1, name="flow")

    # Objective
    m.setObjective(gp.quicksum(var for (i, j), var in flow.items() if i == "s"), sense=-1)

    # Arc-capacity constraints
    m.addConstrs(
        (flow.sum(i, j) <= capacity[i, j] for i, j in arcs), "cap")

    # Flow-conservation constraints
    m.addConstrs(
        (flow.sum(j, '*') == flow.sum('*', j)
            for j in nodes if j not in ('s', 't')), "node")

    # Compute optimal solution
    m.optimize()

    # Print solution
    if m.status == GRB.OPTIMAL:
        solution = m.getAttr('x', flow)
        print('\nOptimal flows')
        for i, j in arcs:
            if solution[i, j] > 0:
                print('%s -> %s: %g' % (i, j, solution[i, j]))
```


```python
max_flow(nodes, arcs, capacity)
```

    Using license file /Library/gurobi/gurobi.lic
    Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
    Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
    Optimize a model with 12 rows, 8 columns and 20 nonzeros
    Model fingerprint: 0xfa101f8c
    Coefficient statistics:
      Matrix range     [1e+00, 1e+00]
      Objective range  [1e+00, 1e+00]
      Bounds range     [0e+00, 0e+00]
      RHS range        [9e+01, 2e+02]
    Presolve removed 12 rows and 8 columns
    Presolve time: 0.01s
    Presolve: All rows and columns removed
    Iteration    Objective       Primal Inf.    Dual Inf.      Time
           0    1.8000000e+02   0.000000e+00   0.000000e+00      0s
    
    Solved in 0 iterations and 0.01 seconds
    Optimal objective  1.800000000e+02
    
    Optimal flows
    s -> A: 90
    s -> B: 90
    A -> C: 90
    B -> D: 90
    C -> t: 90
    D -> t: 90


## Minimum Cut

Next, let's consider the problem of determining the minimum cut.

### Variables

The variables used are as follows

|  Parameter | Type | Description | 
|:---------|:--------| :----- |
| $r_{i,j}$ | Continuous | 1 if arc $(i,j) \in A$ is removed  |
| $z_{i,j}$ | Continuous | 1 if $i \in V \setminus \{s,t\}$ is removed  |

### Model

A model of the problem can then be defined as follows:

$$
\begin{alignat}3
\text{minimize} \ & \sum_{(i,j) \in A} c_{i,j} \cdot r_{i,j} && & (2a) \\
s.t.\ & r_{s,j} + z_j \geq 1 \quad && \forall j \in V : (s, j) \in A & \quad (2c)\\
& r_{i,j} + z_j \geq z_i \quad && \forall (i, j) \in A: i \neq s \text{ and } j \neq t & \quad (2b) \\
& r_{i,t} - z_i \geq 0 \quad && \forall i \in V : (i, t) \in A & \quad (2d)
\end{alignat}
$$

The objective (2a) is to minimise the sum of capacities of the arcs that are removed from the network. Constraints (2b) ensure that for all arcs leaving the source node, either the adjacent node $j$ is connected to the sink, i.e., $z_j = 1$, or the arc is removed, $r_{s, j} = 1$. Constraints (2c) ensure that, for each arc where both nodes are neither the source or the sink, that if predecessor is connected, $z_i = 1$, then either the successor is connected $z_j =1$ or the arc is removed $r_{i,j} = 1$. Finally constraints (2d) ensure that for any arc adjacent to the sink node, if the predecessor is connected, $z_i = 1$, then the arc must be removed, $r_{i,t}=1$. 


```python
def min_cut(nodes, arcs, capacity):

    # Create optimization model
    m = gp.Model('cut')
    m.ModelSense = 1

    # Create variables
    remove = m.addVars(arcs, vtype=GRB.CONTINUOUS, obj=capacity, name="r_")
    connect = m.addVars((i for i in nodes if i not in ("s", "t")), name="z_", vtype=GRB.CONTINUOUS)

    # Arc-capacity constraints
    for (i, j) in arcs:

        if i == "s":
            m.addConstr(remove["s", j] + connect[j] >= 1)
        
        elif j == "t":
            m.addConstr(remove[i, "t"] - connect[i] >= 0)
        
        else:
            m.addConstr(remove[i, j] + connect[j] - connect[i] >= 0)

    # Compute optimal solution
    m.optimize()

    # Print solution
    if m.status == GRB.OPTIMAL:
        solution = m.getAttr('x', remove)
        print('\nOptimal cuts')
        for i, j in arcs:
            if solution[i, j] > 0.5:
                print('%s -> %s: %g' % (i, j, capacity[i, j]))
```


```python
min_cut(nodes, arcs, capacity)
```

    Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
    Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
    Optimize a model with 8 rows, 12 columns and 20 nonzeros
    Model fingerprint: 0x28bf68bd
    Coefficient statistics:
      Matrix range     [1e+00, 1e+00]
      Objective range  [9e+01, 2e+02]
      Bounds range     [0e+00, 0e+00]
      RHS range        [1e+00, 1e+00]
    Presolve removed 6 rows and 9 columns
    Presolve time: 0.01s
    Presolved: 2 rows, 3 columns, 4 nonzeros
    
    Iteration    Objective       Primal Inf.    Dual Inf.      Time
           0    9.0000000e+01   1.000000e+00   0.000000e+00      0s
           1    1.8000000e+02   0.000000e+00   0.000000e+00      0s
    
    Solved in 1 iterations and 0.01 seconds
    Optimal objective  1.800000000e+02
    
    Optimal cuts
    A -> C: 90
    D -> t: 90


### Questions

* How do the number of variables of one model compare with the number of constraints of the other?
* How do the optimal solutions compare?
* How do the objective values of feasible solutions to the max flow problem compare with those of the min cut problem?
* Why are the $r$ and $c$ variables continuous when a cut is clearly a binary operation?

4df2b12ff816be13a4808ab92ec8903849d66096
ipynb
Jupyter Notebook
week2/Max Flow + Min Cut.ipynb
stevedwards/gurobi_course
badb716d5dac86a77712908637cbc08722d0415d
MIT
week2/Max Flow + Min Cut.ipynb
stevedwards/gurobi_course
badb716d5dac86a77712908637cbc08722d0415d
MIT
week2/Max Flow + Min Cut.ipynb
stevedwards/gurobi_course
badb716d5dac86a77712908637cbc08722d0415d
MIT
Qwen/Qwen-72B
1. YES
2. YES