Min Cost Max Flow, in depth: residual costs, the negative cycle condition, successive shortest paths with potentials and modelling

The problem, stated precisely

A network has vertices, directed edges with capacity u(e) and cost c(e) per unit, a source s and a sink t. A flow f assigns 0 ≤ f(e) ≤ u(e) to each edge and conserves flow at every vertex except s and t. Its cost is the sum of f(e) × c(e). Two variants are common and worth separating:

• Min cost flow of value k: send exactly k units from s to t at minimum cost (fails if the maximum flow is below k).
• Min cost max flow: among all maximum flows, find the cheapest. This is the first variant with k set to the max flow value, and the algorithm below solves both by stopping early or not.

A more general form gives each vertex a supply or demand b(v) and asks for the cheapest flow meeting all of them; adding a super source and super sink reduces it to the s-t form. The problem is a linear program whose constraint matrix is totally unimodular, so with integer capacities there is always an integer optimal solution, and the combinatorial algorithms below find one directly. The dual variables of that LP are exactly the potentials used later; LP duality explains the connection.

Residual graphs with costs, and the optimality condition

As in max flow, the residual graph has a forward arc with remaining capacity u(e) − f(e) and a reverse arc with capacity f(e). The new rule is about cost: the forward arc costs c(e), and the reverse arc costs −c(e), because pushing flow back along it cancels flow you previously paid for, and refunds that cost.

This gives the central theorem. A feasible flow of value k has minimum cost if and only if its residual graph contains no negative cost cycle. One direction is immediate: if a negative cycle exists, push one unit around it. Conservation still holds, the flow value is unchanged, and the cost drops, so the flow was not optimal. The other direction follows from decomposing the difference between any two flows of value k into cycles in the residual graph; if none is negative, no other flow is cheaper.

The theorem suggests two algorithms. Cycle cancelling starts from any max flow and repeatedly cancels negative cycles. Successive shortest paths keeps the residual graph free of negative cycles from the start and grows the flow along cheapest paths. The second is simpler and usually faster for moderate flow values.

Successive shortest paths

Successive shortest paths starts from the zero flow. If the original graph has no negative cycles, the zero flow is optimal for value 0. Then it repeats: find a cheapest s-t path in the residual graph, push as much as the path's bottleneck allows, and update residual capacities. It stops when t is unreachable or the target value is reached.

Why does this stay optimal? Suppose the current flow has no negative residual cycle, and let d(v) be shortest path distances from s. Augmenting along a shortest path only creates reverse arcs on edges where d(v) = d(u) + c(u,v), and each reverse arc has reduced cost zero with respect to d, so no negative cycle appears. By induction, every intermediate flow is a minimum cost flow for its value. That is also why the cost per unit of successive paths never decreases: the minimum total cost is a convex function of k.

Potentials: making Dijkstra safe

Residual graphs contain negative arcs (the reverse arcs), so plain Dijkstra is not safe. The fix is the same reweighting Johnson's algorithm uses. Keep a potential π(v) per vertex and define the reduced cost of an arc as c(u,v) + π(u) − π(v). Around any cycle the potentials cancel, so reduced costs preserve which paths are shortest, and if π equals the shortest distances, every residual arc has a non-negative reduced cost.

The loop then becomes: run Dijkstra on reduced costs to get distances d; set π(v) += d(v) for every reachable v; augment. After the update, arcs on the shortest path tree have reduced cost zero and the new reverse arcs also have zero reduced cost, so the invariant survives. The initial potentials come from one Bellman-Ford pass if any original cost is negative, or are all zero otherwise. Bellman-Ford also detects a negative cycle reachable from s, which successive shortest paths cannot handle on its own; a negative cycle elsewhere must be cancelled or rejected before you start.

A tested implementation

The implementation stores edges in adjacency lists with the index of the paired reverse edge, so updating residual capacity is constant time. It returns the flow, the cost, and the (amount, unit cost) of each augmentation for inspection.

import heapq

class MinCostFlow:
    def __init__(self, n):
        self.n = n
        self.g = [[] for _ in range(n)]   # edge = [to, cap, cost, rev_index]

    def add_edge(self, u, v, cap, cost):
        self.g[u].append([v, cap, cost, len(self.g[v])])
        self.g[v].append([u, 0, -cost, len(self.g[u]) - 1])

    def _bellman_ford(self, s):
        INF = float("inf")
        dist = [INF] * self.n
        dist[s] = 0
        for _ in range(self.n - 1):
            changed = False
            for u in range(self.n):
                if dist[u] == INF:
                    continue
                for v, cap, cost, _ in self.g[u]:
                    if cap > 0 and dist[u] + cost < dist[v]:
                        dist[v] = dist[u] + cost
                        changed = True
            if not changed:
                return dist
        for u in range(self.n):
            if dist[u] == INF:
                continue
            for v, cap, cost, _ in self.g[u]:
                if cap > 0 and dist[u] + cost < dist[v]:
                    raise ValueError("negative cycle reachable from source")
        return dist

    def flow(self, s, t, limit=float("inf")):
        INF = float("inf")
        pot = [0 if d == INF else d for d in self._bellman_ford(s)]
        total_flow, total_cost, steps = 0, 0, []
        while total_flow < limit:
            dist = [INF] * self.n
            prev = [None] * self.n            # (node, edge index)
            dist[s] = 0
            pq = [(0, s)]
            while pq:
                d, u = heapq.heappop(pq)
                if d > dist[u]:
                    continue
                for i, (v, cap, cost, _) in enumerate(self.g[u]):
                    if cap <= 0:
                        continue
                    nd = d + cost + pot[u] - pot[v]   # reduced cost, >= 0
                    if nd < dist[v]:
                        dist[v] = nd
                        prev[v] = (u, i)
                        heapq.heappush(pq, (nd, v))
            if dist[t] == INF:
                break
            for v in range(self.n):
                if dist[v] < INF:                     # unreachable nodes keep their potential
                    pot[v] += dist[v]
            push, v = limit - total_flow, t
            while v != s:
                u, i = prev[v]
                push = min(push, self.g[u][i][1])
                v = u
            v = t
            while v != s:
                u, i = prev[v]
                e = self.g[u][i]
                e[1] -= push
                self.g[v][e[3]][1] += push
                v = u
            unit = pot[t] - pot[s]                    # true cost of this path
            steps.append((push, unit))
            total_flow += push
            total_cost += push * unit
        return total_flow, total_cost, steps

This code was checked against exhaustive enumeration of all integer flows on 300 small random graphs without negative cycles, most of them with negative edge costs, for every flow value from zero to the maximum. Do the same for your own implementation: a brute-force oracle on tiny inputs catches the potential bugs listed below faster than any amount of reasoning.

Worked example

Run the algorithm from S to T. Every cost is non-negative, so the initial potentials are zero. Each iteration finds the cheapest residual path:

Now T is unreachable: both edges out of S are saturated, so the maximum flow is 6 and its minimum cost is 49. The unit costs 5, 8, 9, 10, 12 never decrease, as the theory predicts, and the cheapest flow for any smaller target can be read off the table: 1 unit costs 5, 3 units cost 18, 4 units cost 27. In this network no augmentation needed a reverse arc. In tighter networks one often does: a later path cancels part of an earlier one, which is how the algorithm undoes an early greedy choice. Suurballe's algorithm is a two-path special case where that cancellation is the whole point.

Complexity and faster alternatives

Each iteration costs one Dijkstra, O(E log V) with a binary heap, and pushes at least one unit, so with integer capacities and total flow F the algorithm runs in O(F · E log V) after one O(VE) Bellman-Ford. That bound is pseudo-polynomial: it grows with the numeric value F, not just the size of the graph. When F is small, as in assignment problems where F is at most the number of workers, this is excellent. When capacities are large, it can be very slow.

In production, the usual answer is a well-tested library implementation of network simplex or cost scaling. Write your own when the graph is small, when you need to embed it in a hot loop, or when you need incremental behaviour such as the per-unit cost table above.

Modelling problems as min cost flow

Most of the skill is in the modelling. Common patterns:

• Assignment. Source to each worker with capacity 1 and cost 0, worker to job with capacity 1 and the assignment cost, job to sink with capacity 1. Min cost max flow is the optimal assignment, the same problem the Hungarian method solves.
• Transportation. Source to each warehouse with capacity equal to stock, warehouse to store with shipping cost per unit, store to sink with capacity equal to demand.
• Maximising profit. Negate profits to make them costs. The graph now has negative edges, so you need the Bellman-Ford start, and you must decide whether you want maximum flow or maximum profit; for the latter, stop as soon as a path's unit cost becomes non-negative.
• Convex costs. If the cost of sending x units is convex and piecewise linear, split the edge into parallel edges with increasing unit costs; the algorithm fills the cheap ones first. Concave costs, such as economies of scale, cannot be expressed this way and make the problem NP-hard in general.
• Lower bounds and demands. Edges that must carry at least some flow reduce to demands at their endpoints, then to an s-t problem with a super source and sink.

Implementation bugs and failure modes

• Reverse arc with the wrong cost. Giving reverse arcs cost +c instead of −c produces plausible but non-optimal answers. A brute-force oracle exposes it instantly.
• Updating unreachable potentials. Adding an infinite distance to a potential poisons it, and later reduced costs become infinite or NaN. Update only vertices Dijkstra reached.
• Dijkstra without potentials. Works on the first iteration, then silently returns wrong paths once reverse arcs with negative cost appear.
• Negative cycles in the input. Successive shortest paths assumes the zero flow is optimal. If the input has a negative cycle, detect it and either cancel it first or report an error.
• Floating-point costs. Reduced costs that should be zero drift slightly negative and Dijkstra's assumptions break. Scale costs to integers where possible.
• Confusing the two objectives. Code that stops at a flow limit returns the cheapest flow of that value, not the cheapest maximum flow. State which one a caller gets.

What to do next

1. Implement the class above and reproduce the worked example: flow 6, cost 49, unit costs 5, 8, 9, 10, 12.
2. Write a brute-force oracle that enumerates integer flows on graphs with up to six edges and compare on random inputs, including negative costs.
3. Model a small assignment problem with five workers and five jobs and check it against a Hungarian implementation.
4. Add a convex cost edge as parallel arcs and confirm the cheap arcs fill first.
5. Time your implementation as total flow grows, to see the pseudo-polynomial bound in practice.
6. Read about network simplex and cost scaling, then compare your code with a library solver on a graph with large capacities.
7. Review Dinic's algorithm for the max flow side, since many min cost flow problems start as a max flow model.