Maximum Flow, in depth: modelling, solver libraries, reading the minimum cut and re-solving when capacity changes

The model from first principles

A flow network is a directed graph with a source s, a sink t and a non-negative capacity c(u, v) on every arc. A flow assigns each arc a value f(u, v) obeying 0 <= f(u, v) <= c(u, v), and every vertex except the source and sink passes on exactly what it receives. Maximum flow asks for the largest net amount that can leave the source.

Two ideas make the problem tractable. The first is the residual graph: each arc carrying flow gets a reverse arc whose capacity is the flow sent, so a later step can undo a greedy choice. A flow is maximum exactly when the residual graph has no path from source to sink. The second is the cut: split the vertices into a set containing the source and a set containing the sink, and add up the capacities of arcs crossing from the first to the second. No flow can exceed any cut, and the max-flow min-cut theorem says the largest flow equals the smallest cut. That equality is the practical payoff. The maximum flow tells you how much you can do, and the minimum cut tells you why you cannot do more.

A third fact matters for modelling: with integer capacities, the standard algorithms return an integer maximum flow, which is what lets flow answer discrete questions such as which job goes to which rack.

The algorithm landscape

You rarely need to choose an algorithm from scratch, but you should know what the library is doing. The families below all compute the same answer; they differ in how they find it.

Worst-case bounds are loose on real inputs. Push-relabel with global-relabel and gap heuristics is the usual library workhorse: networkx defaults to preflow_push, and OR-Tools ships a push-relabel solver behind SimpleMaxFlow.

Worked example: placing GPU jobs on racks

A scheduler has four jobs waiting. J1 and J4 need H100s, J3 needs A100s and J2 runs on either. J1 wants 10 GPUs, J2 wants 6, J3 and J4 want 4 each: demand 24. H100 racks R1 and R2 have 8 and 4 free GPUs, A100 rack R3 has 12. Supply is also 24, so a naive check says everything fits.

Build the network in three layers. An arc from the source to each job has capacity equal to its demand. An arc from a job to each rack it is allowed to use carries the job's demand as its capacity. An arc from each rack to the sink has the rack's free GPUs as its capacity. One unit of flow is one GPU placed, so the maximum flow is the largest number of GPUs that can be placed while respecting both demand and eligibility.

The solver returns 22, not 24. The minimum cut explains the gap without any further analysis: the cut consists of both H100 rack arcs (8 + 4) and the source arcs of the two jobs that can use A100s (6 + 4). In words, the H100-only jobs need 14 GPUs and only 12 H100s exist. The A100 rack keeps 2 GPUs idle. Adding A100 capacity cannot help; freeing two H100s, or relaxing J1 or J4 to run on A100, would.

Solving it with a library

Do not write your own solver for a one-off analysis. In networkx, nx.maximum_flow returns both the value and a nested dict of per-arc flows.

import networkx as nx

JOBS = {"J1": (10, {"h100"}), "J2": (6, {"h100", "a100"}),
        "J3": (4, {"a100"}), "J4": (4, {"h100"})}
RACKS = {"R1": ("h100", 8), "R2": ("h100", 4), "R3": ("a100", 12)}

def build(jobs, racks):
    g = nx.DiGraph()
    for j, (need, kinds) in jobs.items():
        g.add_edge("s", j, capacity=need)
        for r, (kind, free) in racks.items():
            if kind in kinds:
                g.add_edge(j, r, capacity=need)
    for r, (kind, free) in racks.items():
        g.add_edge(r, "t", capacity=free)
    return g

g = build(JOBS, RACKS)
value, flow = nx.maximum_flow(g, "s", "t")          # 22
cut_value, (S, T) = nx.minimum_cut(g, "s", "t")     # 22
placement = {j: {r: f for r, f in flow[j].items() if f} for j in JOBS}

For a service that solves repeatedly, a compiled solver is worth the extra ceremony. OR-Tools numbers vertices as integers and returns arc indices you keep for reading flows back:

from ortools.graph.python import max_flow

idx = {name: i for i, name in enumerate(g.nodes)}
smf = max_flow.SimpleMaxFlow()
arcs = [(u, v, smf.add_arc_with_capacity(idx[u], idx[v], d["capacity"]))
        for u, v, d in g.edges(data=True)]
if smf.solve(idx["s"], idx["t"]) != smf.OPTIMAL:
    raise RuntimeError("max flow did not solve")
print(smf.optimal_flow())                            # 22
names = {i: n for n, i in idx.items()}
print(sorted(names[i] for i in smf.get_source_side_min_cut()))
# ['J1', 'J4', 'R1', 'R2', 's']

Run both and you will notice that networkx returned a different minimum cut, with source side {s, J1, J2, J3, J4, R1, R2}, cutting R1, R2 and the two job-to-R3 arcs. It also totals 22. Minimum cuts are often not unique. The set of vertices reachable from the source in the final residual graph is the smallest source side and the most useful for explanations, so compute it yourself rather than trusting whichever cut a library reports.

Reading the answer: value, flows and cut

A flow answer has three parts, and production code should use all of them.

• The value is the headline: 22 of 24 GPUs placed.
• The per-arc flows are the decision: J1 got 4 GPUs on R1 and 4 on R2, J4 got 4 on R1, J2 and J3 got 6 and 4 on R3.
• The cut is the explanation. Each arc in the minimal cut is a binding constraint; translated to domain terms, H100 capacity is short by two.

Notice what the flow did not do. J1 asked for 10 GPUs and received 8, split across two racks, but a real training job usually needs all of its GPUs at once. Max flow cannot express all-or-nothing constraints; adding them makes the problem an integer program. Use flow as a fast feasibility and diagnosis step, then settle all-or-nothing decisions with a heuristic or an exact solver.

Before acting on any answer, verify it. Check capacity and conservation on every arc and vertex, then confirm that the flow value equals the capacity of the cut you computed. Equal values prove optimality in linear time, which is cheap insurance against a modelling bug or a solver called on the wrong graph. The min-cut theorem article walks through why that certificate is sufficient.

Re-solving when capacity changes

Clusters change constantly: a rack is drained, a link degrades, a job is cancelled. For large graphs the previous flow is a valuable starting point, and most library interfaces discard it. Keeping your own residual graph lets you repair instead of restart.

Increasing a capacity never invalidates the current flow, so you simply keep augmenting. Decreasing a capacity below the flow on that arc breaks the capacity rule, and the repair takes three steps: remove the excess from the arc, try to reroute it around the arc in the residual graph, and if that fails, cancel the remainder by pushing it back towards the source on one side and from the sink on the other. Then augment from source to sink as usual.

from collections import deque

class IncrementalMaxFlow:
    """Edmonds-Karp on arc arrays that keep their flow between solves."""

    def __init__(self, n):
        self.n, self.adj = n, [[] for _ in range(n)]
        self.to, self.cap, self.flow = [], [], []

    def add_edge(self, u, v, c):
        for a, b, cc in ((u, v, c), (v, u, 0)):
            self.adj[a].append(len(self.to))
            self.to.append(b); self.cap.append(cc); self.flow.append(0)
        return len(self.to) - 2        # forward arc id; its twin is id ^ 1

    def augment(self, s, t, limit=float("inf")):
        pushed = 0
        while s != t and pushed < limit:
            prev, q = {s: None}, deque([s])
            while q and t not in prev:
                u = q.popleft()
                for e in self.adj[u]:
                    v = self.to[e]
                    if v not in prev and self.cap[e] > self.flow[e]:
                        prev[v] = e
                        q.append(v)
            if t not in prev:
                break
            path, v = [], t
            while v != s:
                path.append(prev[v]); v = self.to[prev[v] ^ 1]
            d = min(min(self.cap[e] - self.flow[e] for e in path), limit - pushed)
            for e in path:
                self.flow[e] += d; self.flow[e ^ 1] -= d
            pushed += d
        return pushed

    def value(self, s):
        return sum(self.flow[e] for e in self.adj[s])

    def set_capacity(self, e, new_cap, s, t):
        self.cap[e] = new_cap
        excess = self.flow[e] - new_cap
        if excess > 0:
            u, v = self.to[e ^ 1], self.to[e]
            self.flow[e] -= excess; self.flow[e ^ 1] += excess
            left = excess - self.augment(u, v, excess)    # reroute around the arc
            if left:
                self.augment(u, s, left)                   # return inflow to the source
                self.augment(t, v, left)                   # and pull outflow back from the sink
        self.augment(s, t)
        return self.value(s)

On the worked example, draining R1 from 8 free GPUs to 2 drops the answer from 22 to 16, and then shrinking R3 from 12 to 9 drops it to 15; both match a fresh networkx solve. Test a repair routine the way this one was tested: 400 random graphs, five capacity changes each, comparing against a from-scratch solve and asserting capacity and conservation after every change. A missed cancellation leaves a flow that looks plausible and is infeasible.

Operational guidance

Use integers. Express capacities in whole units: GPUs, megabits per second, minutes. With floats a residual of 1e-12 looks like spare capacity and solvers can crawl through near-zero augmentations. If you must use reals, scale and round.

Watch the graph size. A thousand jobs each eligible for a thousand machines is a million arcs. Group interchangeable machines into one vertex per pool, then expand afterwards.

Know when flow is the wrong tool. If arcs have costs and you want the cheapest way to place everything, you need minimum-cost flow. If some arcs need a minimum amount, see flows with lower bounds and demands. If decisions are all-or-nothing, you are in integer programming territory.

Failure modes

• Undirected edges modelled as one arc. A cable carries traffic both ways; model it as two opposite arcs with the full capacity each, or the answer will be too small.
• Vertex capacities ignored. If a switch or a reviewer has a limit, split the vertex into an in-vertex and an out-vertex joined by an arc with that capacity.
• Trusting a solver's cut for explanations. As the worked example showed, the cut returned may not be the minimal one; compute reachability from the source yourself.

Trade-offs

What to do next

1. Pick one allocation problem you own and draw it as source, layers and sink, writing the unit of one unit of flow.
2. Solve it with networkx and print the value, the per-arc flows and the source side of the cut.
3. Translate each saturated cut arc into a sentence a human would act on, and check it against intuition.
4. Add a verifier that checks capacity, conservation and flow equals cut capacity before any decision is used.
5. If the graph changes often, prototype an incremental solver and randomise-test it against full solves.
6. Read Dinic or push-relabel before tuning performance.