vii. Min-Cost Flow

Edges have capacities and per-unit costs; send a required amount (or the maximum) from $s$ to $t$ as cheaply as possible, or more generally route given supplies and demands with lower bounds. A flow is optimal iff its residual graph has no negative cycle, iff there are potentials $\pi$ with reduced costs $c(u,v)+\pi (u)-\pi (v)\ge 0$ on every residual edge.

Method	Time	Use
successive shortest paths (SSP) + potentials	$O(F\cdot E\log V)$	the default; $F$ is usually small (matchings)
SSP with $O(V^2)$ Dijkstra	$O(F{V}^2)$	dense graphs
capacity scaling ($b$-flow)	$O(E^2\log U\log V)$	huge capacities, supplies, negative cycles
cost scaling (push–relabel)	$O(V^{2}E\log (VC))$	large graphs, any costs
network simplex	no good bound, fastest in practice	large graphs, any costs
min-mean cycle cancelling	strongly polynomial	theory

Successive shortest paths

Augment along a cheapest residual path each time. Dijkstra works on reduced costs once potentials are valid: start from Bellman–Ford if some cost is negative (no negative cycle allowed), else from zero, and add the distances to the potentials after each round.

struct MCMF {
    struct E { int to; ll cap, cost; };
    int n; vector<E> e; vector<vector<int>> g; vector<ll> pot, dist; vector<int> pe;
    MCMF(int n) : n(n), g(n), pot(n), dist(n), pe(n) {}
    int addEdge(int u, int v, ll cap, ll cost) {
        g[u].push_back(sz(e)); e.push_back({v, cap, cost});
        g[v].push_back(sz(e)); e.push_back({u, 0, -cost});
        return sz(e) - 2;
    }
    void initPotentials(int s) {                               // only needed with negative costs
        const ll INF = LLONG_MAX; fill(all(pot), INF); pot[s] = 0;
        for (int it = 0; it < n; it++) {
            bool ch = false;
            for (int u = 0; u < n; u++) if (pot[u] != INF)
                for (int id : g[u]) if (e[id].cap > 0 && pot[u] + e[id].cost < pot[e[id].to]) pot[e[id].to] = pot[u] + e[id].cost, ch = true;
            if (!ch) break;
        }
        for (ll& p : pot) if (p == INF) p = 0;                 // unreachable vertices stay unreachable
    }
    pair<ll, ll> run(int s, int t, ll need = LLONG_MAX) {      // {flow, cost}
        ll flow = 0, cost = 0; const ll INF = LLONG_MAX;
        while (flow < need) {
            fill(all(dist), INF); dist[s] = 0;
            priority_queue<pair<ll, int>, vector<pair<ll, int>>, greater<>> pq; pq.push({0, s});
            while (sz(pq)) {
                auto [d, u] = pq.top(); pq.pop();
                if (d > dist[u]) continue;
                for (int id : g[u]) {
                    if (e[id].cap <= 0) continue;
                    ll nd = d + e[id].cost + pot[u] - pot[e[id].to];
                    if (nd < dist[e[id].to]) dist[e[id].to] = nd, pe[e[id].to] = id, pq.push({nd, e[id].to});
                }
            }
            if (dist[t] == INF) break;
            for (int v = 0; v < n; v++) if (dist[v] < INF) pot[v] += dist[v];
            ll push = need - flow;
            for (int v = t; v != s; v = e[pe[v] ^ 1].to) push = min(push, e[pe[v]].cap);
            for (int v = t; v != s; v = e[pe[v] ^ 1].to) e[pe[v]].cap -= push, e[pe[v] ^ 1].cap += push, cost += push * e[pe[v]].cost;
            flow += push;
        }
        return {flow, cost};
    }
};

• Cheapest flow of any value: stop as soon as the shortest path has non-negative cost ($\text{pot}[t]-\text{pot}[s]\ge 0$ after the update).
• The cost of the first $k$ units is convex in $k$ (path costs are non-decreasing), so SSP answers “minimum cost for every $k$” in one run, and it plugs into Aliens / slope tricks.
• Max-cost flow: negate costs (then initialise potentials).

Proof: optimality criterion and why SSP is correct

Negative cycle ⟹ not optimal: pushing $\epsilon$ around a negative residual cycle keeps conservation and the flow value and lowers the cost. No negative cycle ⟹ optimal: if $f^{\ast }$ is cheaper with the same value, $f^{\ast }-f$ is a circulation in the residual graph of $f$, decomposing into residual cycles; their total cost is negative, so one is negative. Potentials: with no negative cycle, shortest distances $d$ exist and satisfy $d(v)\le d(u)+c(u,v)$ on residual edges, i.e. reduced costs w.r.t. $d$ are $\ge 0$; conversely non-negative reduced costs make every cycle’s cost (which equals its reduced cost) non-negative.

SSP keeps the invariant: suppose reduced costs w.r.t. $\pi$ are $\ge 0$ and $d$ are Dijkstra distances under them. With ${\pi }^{\prime }=\pi +d$, every residual edge has $c^{{\pi }^{\prime }}(u,v)=c^{\pi }(u,v)+d(u)-d(v)\ge 0$, and edges on the shortest path have $c^{{\pi }^{\prime }}=0$. Augmenting creates reverse edges of those path edges, with reduced cost $0$. So reduced costs stay non-negative, the new residual graph has no negative cycle, and by induction each intermediate flow is a cheapest flow of its value. Convexity: distances under the invariant never decrease between rounds (the new reverse edges have reduced cost 0 and lie on old shortest paths), so the cost per unit $\text{pot}[t]-\text{pot}[s]$ is non-decreasing. Termination with integral capacities: each round adds $\ge 1$ unit.

Modelling with costs

Requirement	Construction
assignment / weighted bipartite matching	Assignment Problem (Hungarian is faster)
transportation (supplies, demands, shipping costs)	$s\to$ supplier (supply, 0), supplier $\to$ consumer (cap, cost), consumer $\to t$ (demand, 0)
convex cost $g(x)$ on an edge	unit-capacity parallel edges with costs $g(1)-g(0),g(2)-g(1),\dots$ (increasing, so they fill in order)
$k$ disjoint paths of minimum total length	unit capacities (split vertices for vertex-disjoint), flow $k$ (Suurballe)
interval scheduling with weights, $k$ machines	chain $0\to 1\to \cdots \to T$ with capacity $k$ and cost 0, each interval $l\to r$ with capacity 1 and cost $-w$
reward for a choice	negative cost; initialise potentials (or shift costs by a constant per unit if every path uses one such edge)
lower bounds, supplies, negative cycles	min-cost $b$-flow below

Proof: convex costs via parallel edges

Split an edge with convex $g$ into unit edges of costs ${\delta }_i=g(i)-g(i-1)$, which are non-decreasing by convexity. For a fixed total $x$ through the bundle, the cheapest way to choose $x$ of these unit edges is the $x$ cheapest, costing ${\sum }_{i\le x}{\delta }_i=g(x)-g(0)$, and an optimal flow always picks the cheapest ones (otherwise swap a used expensive unit for an unused cheaper one). So the optimum over the split network equals the optimum with the convex edge costs. With non-convex $g$ the cheapest units would be the “later” ones, which is not a valid representation.

Min-cost $b$-flow (capacity scaling, negative cycles allowed)

Supplies $b_v$ (positive = produce), edges with $\ell \le f\le u$ and any integer costs. Each vertex is joined to a virtual root by an edge with cost $-\infty$ (big-M) to force supplies through; capacities are revealed bit by bit, and each new $\Delta$ is pushed along shortest paths with potentials.

struct MinCostBFlow {
    struct E { int src, dst; ll flow, cap, cost; };
    int N; ll INF = 1; i128 base = 0;
    vector<ll> B, P, lower; vector<vector<int>> G; vector<E> e; vector<int> pre;
    MinCostBFlow(int n) : N(n), B(n), P(n + 1), G(n + 1), pre(n + 1) {}
    int ae(int a, int b, ll l, ll u, ll c) {
        G[a].push_back(sz(e)); e.push_back({a, b, 0, u - l, c});
        G[b].push_back(sz(e)); e.push_back({b, a, 0, 0, -c});
        base += (i128)l * c;
        if (a < N) B[a] -= l;
        if (b < N) B[b] += l;
        return sz(e) - 2;
    }
    int addEdge(int u, int v, ll lo, ll hi, ll cost) { lower.push_back(lo); return ae(u, v, lo, hi, cost) / 2; }
    void addSupply(int v, ll b) { B[v] += b; }
    void pin(int x) {                                         // shortest paths from x, folded into P
        priority_queue<pair<ll, int>, vector<pair<ll, int>>, greater<>> pq; vector<ll> add(N + 1);
        for (int i = 0; i <= N; i++) add[i] = -P[i];
        add[x] -= 3 * INF; pre[x] = -1;
        for (int i = 0; i <= N; i++) pq.push({add[i], i});
        while (sz(pq)) {
            auto [d, u] = pq.top(); pq.pop();
            if (d > add[u]) continue;
            for (int id : G[u]) if (e[id].flow < e[id].cap) {
                ll nd = d + e[id].cost + P[e[id].src] - P[e[id].dst];
                if (nd < add[e[id].dst]) add[e[id].dst] = nd, pre[e[id].dst] = id, pq.push({nd, e[id].dst});
            }
        }
        for (int i = 0; i <= N; i++) P[i] += add[i];
    }
    void addCap(int id, ll d) {
        if (e[id].flow == e[id].cap) {
            pin(e[id].dst);
            if (P[e[id].src] - P[e[id].dst] + e[id].cost < 0) {   // closing a negative cycle: push d around it
                for (int x = e[id].src; pre[x] != -1; x = e[pre[x]].src) e[pre[x]].flow += d, e[pre[x] ^ 1].flow -= d;
                e[id].flow += d, e[id ^ 1].flow -= d;
            }
        }
        e[id].cap += d;
    }
    optional<i128> solve() {                                  // nullopt if infeasible
        int M = sz(e);
        for (int i = 0; i < M; i += 2) INF += abs(e[i].cost);
        for (int i = 0; i < N; i++) { int a = N, b = i; ll c = B[i]; if (c < 0) c = -c, swap(a, b); ae(a, b, 0, c, -INF); }
        vector<ll> caps(sz(e)); ll mx = 1;
        for (int i = 0; i < sz(e); i++) swap(caps[i], e[i].cap), mx = max(mx, caps[i]);
        ll d = 1; while (d < mx) d *= 2;
        for (; d; d /= 2) for (int i = 0; i < sz(e); i++) if (caps[i] & d) addCap(i, d);
        for (int i = 0; i < N; i++) if (e[M + 2 * i].flow < e[M + 2 * i].cap) return nullopt;
        i128 ans = base;
        for (int i = 0; i < M; i += 2) ans += (i128)e[i].flow * e[i].cost;
        return ans;
    }
    ll flowOn(int id) const { return lower[id] + e[2 * id].flow; }
    ll potential(int v) const { return P[v]; }
};

$s$–$t$ min-cost flow with negative cycles allowed elsewhere in the graph: add an edge $t\to s$ with cost $-M$ (larger in magnitude than any simple path cost) and capacity up to the max flow value, turning the problem into a pure circulation (zero supplies/demands) that this structure solves directly. (Port of Benq’s capacity scaling, verified on Library Checker “Minimum Cost b-flow”.)

Proof sketch: capacity scaling keeps optimality

Invariant: the current flow (with the currently revealed capacities) has no negative residual cycle, certified by potentials $P$. Raising the capacity of edge $(a,b)$ by $d$ only adds residual capacity on that edge. If it was not saturated, the edge already had reduced cost $\ge 0$ and nothing changes. Otherwise pin(b) recomputes potentials as shortest distances from $b$ (all residual reduced costs non-negative again; the $-3\cdot$INF offset makes $b$ the source of the shortest-path tree). If the reduced cost of $(a,b)$ is negative, the cheapest residual path $b⇝a$ plus the edge is the most negative cycle through it, and pushing $d$ around it (all its residual capacities are multiples of $d$ by the scaling invariant) removes it, restoring the invariant. At the end all real capacities are revealed. The big-M edges to the root cost $-\text{INF}$ with INF larger than any simple path cost, so an optimal flow saturates all of them whenever a feasible routing exists, and unsaturated ones certify infeasibility.

Network simplex

Keep a spanning tree basis (rooted at a virtual vertex joined to all by big-M edges); repeatedly bring in an edge with negative reduced cost (block pivoting over $V/3$ edges), push around the cycle it closes in the tree, and drop the bottleneck edge. The fastest practical min-cost flow; same interface as above.

struct NetworkSimplex {
    struct E { int src, dst; ll flow, cap, cost; };
    ll INF = 1; i128 delta = 0; int n;
    vector<ll> B, P, lower; vector<E> e; vector<int> pei, depth; vector<set<int>> tree;
    NetworkSimplex(int n) : n(n), B(n + 1), P(n + 1), pei(n + 1, -1), depth(n + 1), tree(n + 1) {}
    int ae(int a, int b, ll l, ll u, ll c) {
        e.push_back({a, b, 0, u - l, c}); e.push_back({b, a, 0, 0, -c});
        delta += (i128)l * c; B[b] += l, B[a] -= l;
        return sz(e) - 2;
    }
    int addEdge(int u, int v, ll lo, ll hi, ll cost) { lower.push_back(lo); return ae(u, v, lo, hi, cost) / 2; }
    void addSupply(int v, ll b) { B[v] += b; }
    void upd(int ei) { P[e[ei].dst] = P[e[ei].src] + e[ei].cost; pei[e[ei].dst] = ei ^ 1; depth[e[ei].dst] = 1 + depth[e[ei].src]; dfs(e[ei].dst); }
    void dfs(int v) { for (int ei : tree[v]) if (ei != pei[v]) upd(ei); }
    template <class F> void walk(int ei, F f) {               // ei, then the tree path closing its cycle
        f(ei);
        for (int a = e[ei].src, b = e[ei].dst; a != b;) {
            if (depth[a] > depth[b]) f(pei[a] ^ 1), a = e[pei[a]].dst;
            else f(pei[b]), b = e[pei[b]].dst;
        }
    }
    optional<i128> solve() {
        int m = sz(e);
        for (int i = 0; i < m; i += 2) INF += abs(e[i].cost);
        for (int i = 0; i < n; i++) {
            int a = n, b = i; ll c = B[i]; if (c < 0) c = -c, swap(a, b);
            int ei = ae(a, b, 0, c, -INF); tree[a].insert(ei), tree[b].insert(ei ^ 1);
        }
        dfs(n);
        i128 ans = delta; int ptr = 0; const int BLOCK = n / 3 + 1;
        for (int z = 0; z < sz(e) / BLOCK + 1; z++) {
            pair<ll, int> pin{0, -1};
            for (int t = 0; t < BLOCK; ++t, (++ptr) %= sz(e))
                if (e[ptr].flow < e[ptr].cap) pin = min(pin, make_pair(e[ptr].cost + P[e[ptr].src] - P[e[ptr].dst], ptr));
            auto [cost, ein] = pin;
            if (cost == 0) continue;
            pair<ll, int> pout{e[ein].cap - e[ein].flow, ein};
            walk(ein, [&](int ei) { pout = min(pout, make_pair(e[ei].cap - e[ei].flow, ei)); });
            auto [fl, eout] = pout;
            walk(ein, [&](int ei) { e[ei].flow += fl, e[ei ^ 1].flow -= fl; });
            tree[e[ein].src].insert(ein), tree[e[ein].dst].insert(ein ^ 1);
            tree[e[eout].src].erase(eout), tree[e[eout].dst].erase(eout ^ 1);
            upd(pei[e[eout].src] == eout ? ein : ein ^ 1);
            ans += (i128)fl * cost; z = -1;
        }
        for (int i = 0; i < n; i++) {
            if (e[m + 2 * i].flow < e[m + 2 * i].cap) return nullopt;
            ans += (i128)e[m + 2 * i].flow * INF;
        }
        return ans;
    }
    ll flowOn(int id) const { return lower[id] + e[2 * id].flow; }
    ll potential(int v) const { return P[v]; }
};

Proof sketch: network simplex

A basis is a spanning tree $T$ such that every non-tree edge is at its lower bound (flow 0 in the residual sense) or saturated; the tree flow is then determined by conservation. Potentials are defined so that tree edges have reduced cost 0 (upd sets $P(\text{dst})=P(\text{src})+c$ along the tree). If no non-tree residual edge has negative reduced cost, all residual reduced costs are $\ge 0$ and the flow is optimal by the potential criterion. Otherwise the entering edge plus the tree path forms a cycle whose cost equals the entering edge’s reduced cost $<0$ (tree edges contribute 0); pushing the bottleneck amount keeps feasibility, does not increase the cost, and the bottleneck edge leaves, giving a new spanning tree. The big-M root edges make the initial tree feasible. Degenerate pivots (amount 0) can cycle in theory; block pivoting with a rotating pointer avoids it in practice.

Cost scaling

Goldberg–Tarjan: costs are multiplied by $n$, and an $\epsilon$-optimal flow (all residual reduced costs $\ge -\epsilon$) is refined while $\epsilon$ is divided by 4; $\epsilon <1$ after scaling means exactly optimal. First a push–relabel max flow, then the refinement optimises the circulation part too, so negative cycles are fine. Keep $n\cdot \max |c|\cdot F$ well inside 64 bits.

struct CostScaling {                                          // min-cost max-flow, any integer costs
    struct E { int to, rev; ll c, f; };                       // f = residual capacity, c = cost * n
    int N, S, T; ll eps = 0; vector<vector<E>> G;
    vector<unsigned> isq, cur; vector<ll> ex, h; vector<vector<int>> hs; vector<int> co;
    CostScaling(int n, int s, int t) : N(n), S(s), T(t), G(n) {}
    void addEdge(int a, int b, ll cap, ll cost) {
        if (a == b) return;                                   // a non-negative self-loop never helps
        cost *= N; eps = max(eps, abs(cost));
        G[a].push_back({b, sz(G[b]), cost, cap}); G[b].push_back({a, sz(G[a]) - 1, -cost, 0});
    }
    void addFlow(E& e, ll f) {
        E& bk = G[e.to][e.rev];
        if (!ex[e.to] && f) hs[h[e.to]].push_back(e.to);
        e.f -= f; ex[e.to] += f; bk.f += f; ex[bk.to] -= f;
    }
    ll maxFlow() {                                            // highest-label push-relabel
        ex.assign(N, 0); h.assign(N, 0); hs.assign(2 * N, {}); co.assign(2 * N, 0); cur.assign(N, 0);
        h[S] = N; ex[T] = 1; co[0] = N - 1;
        for (auto& e : G[S]) addFlow(e, e.f);
        if (sz(hs[0])) for (int hi = 0; hi >= 0;) {
            int u = hs[hi].back(); hs[hi].pop_back();
            while (ex[u] > 0) {
                if (cur[u] == G[u].size()) {
                    h[u] = 1e9;
                    for (unsigned i = 0; i < G[u].size(); i++) { auto& e = G[u][i]; if (e.f && h[u] > h[e.to] + 1) h[u] = h[e.to] + 1, cur[u] = i; }
                    if (++co[h[u]], !--co[hi] && hi < N) for (int i = 0; i < N; i++) if (hi < h[i] && h[i] < N) --co[h[i]], h[i] = N + 1;
                    hi = h[u];
                } else if (G[u][cur[u]].f && h[u] == h[G[u][cur[u]].to] + 1) addFlow(G[u][cur[u]], min(ex[u], G[u][cur[u]].f));
                else ++cur[u];
            }
            while (hi >= 0 && hs[hi].empty()) --hi;
        }
        return -ex[S];
    }
    void push(E& e, ll amt) { amt = min(amt, e.f); e.f -= amt; ex[e.to] += amt; G[e.to][e.rev].f += amt; ex[G[e.to][e.rev].to] -= amt; }
    void relabel(int v) {
        ll nh = LLONG_MIN / 2;
        for (unsigned i = 0; i < G[v].size(); i++) { const E& e = G[v][i]; if (e.f && nh < h[e.to] - e.c) nh = h[e.to] - e.c, cur[v] = i; }
        h[v] = nh - eps;
    }
    pair<ll, ll> minCostMaxFlow() {
        ll cost = 0;
        for (int i = 0; i < N; i++) for (E& e : G[i]) cost += e.c * e.f;
        ll flow = maxFlow();
        h.assign(N, 0); ex.assign(N, 0); isq.assign(N, 0); cur.assign(N, 0); queue<int> q;
        for (; eps; eps >>= 2) {
            fill(all(cur), 0);
            for (int i = 0; i < N; i++) for (auto& e : G[i]) if (h[i] + e.c - h[e.to] < 0 && e.f) push(e, e.f);
            for (int i = 0; i < N; i++) if (ex[i] > 0) q.push(i), isq[i] = 1;
            while (sz(q)) {
                int u = q.front(); q.pop(); isq[u] = 0;
                while (ex[u] > 0) {
                    if (cur[u] == G[u].size()) relabel(u);
                    for (unsigned& i = cur[u], mx = G[u].size(); i < mx; i++) {
                        E& e = G[u][i];
                        if (h[u] + e.c - h[e.to] < 0) {
                            push(e, ex[u]);
                            if (ex[e.to] > 0 && !isq[e.to]) q.push(e.to), isq[e.to] = 1;
                            if (!ex[u]) break;
                        }
                    }
                }
            }
            if (eps > 1 && eps >> 2 == 0) eps = 1 << 2;
        }
        for (int i = 0; i < N; i++) for (E& e : G[i]) cost -= e.c * e.f;
        return {flow, cost / 2 / N};
    }
    ll flowOn(const E& e) const { return G[e.to][e.rev].f; }
};

Proof sketch: cost scaling

A circulation is $\epsilon$-optimal if some potentials make every residual reduced cost $\ge -\epsilon$. Any circulation is $\max |c|$-optimal (potentials 0). Refinement saturates every residual edge of negative reduced cost (making the pseudo-flow $0$-optimal but unbalanced), then pushes excess only along admissible edges (reduced cost $<0$) and relabels a vertex by lowering its potential to $\max (h(w)-c)-\epsilon$; pushes create reverse edges with positive reduced cost and relabels keep every residual edge at $\ge -\epsilon$, so the result is an $\epsilon$-optimal circulation, with $O(n)$ relabels per vertex by the standard push–relabel argument. With costs multiplied by $n$, a cycle in an $\epsilon$-optimal circulation with $\epsilon <1$ has (scaled) cost $>-n$; scaled costs are multiples of $n$, so it is $\ge 0$: no negative cycle, hence optimal. The max flow computed first fixes the flow value, and refinement only adds circulations, so the result is a minimum-cost maximum flow.