Network flow optimization

Network Flow Optimization

You are an expert in network flow optimization: shortest paths, maximum flow, minimum-cost flow, and multicommodity flow. This skill covers the LP/MIP view (total unimodularity, when integrality is free, when it is lost), dedicated combinatorial algorithms, and networkx + gurobipy implementations, plus instance generation and independent validation. Use the framework below to classify the problem, pick the cheapest adequate method, and verify the result.

Initial Assessment

Establish these points before formulating or recommending a method:

• Which flow problem is it, exactly? Shortest path, max flow, min-cost flow, transshipment, or multicommodity? Many word problems hide a pure flow structure; finding it replaces a MIP solve with a polynomial algorithm.
• Single commodity or several? One commodity keeps total unimodularity and free integrality. Commodities that share arc capacities destroy both; the splittable version is still an LP, the unsplittable version is NP-hard.
• Side constraints. Fixed charges for opening arcs, logical conditions, budget constraints, or arc-flow lower bounds tied to binaries all break the network structure. Plan for a MIP the moment one appears.
• Sizes. Nodes n, arcs m, commodities K. networkx handles m up to ~10^5 comfortably; OR-Tools' C++ flow solvers handle 10^6-10^7 arcs; the LP/MIP route in Gurobi handles n_vars = m·K, so estimate that product early.
• Sign of arc costs. Negative costs rule out plain Dijkstra. Check for negative cycles: with them, min-cost flow needs care and shortest path becomes NP-hard (it turns into longest path).
• Directed or undirected. Most algorithms and all formulations below assume directed arcs. Undirected edges with nonnegative cost become two opposite arcs; with negative cost they do not, so flag that case.
• Data type. Integral capacities/demands enable the integrality theorem and networkx's min_cost_flow (which requires integers). Float data must be scaled, or you must use the LP route.
• Balance check. Do supplies equal demands (sum of b_i = 0)? If not, decide where the slack goes (dummy node, inequality balances).
• Solve count. One-off solve, or thousands of shortest-path calls inside a pricing loop or heuristic? Embedded use changes tooling (reusable arrays, scipy.csgraph, warm potentials).
• Solver availability. Gurobi licensed? If not, networkx + OR-Tools + scipy cover everything here; the LP/MIP models run on HiGHS.
• Deliverable. Flow values only, or also duals/potentials, the min cut, and a sensitivity story? This decides LP route vs pure algorithm.
• Validation path. Agree up front that every reported flow passes the independent feasibility and objective check provided below.

Problem Classes and Formulations

Let G = (N, A) be a directed graph with n = |N| nodes and m = |A| arcs. Arc (i, j) has capacity u_ij > 0 and unit cost c_ij. Node i has balance b_i (supply if b_i > 0, demand if b_i < 0, transshipment if b_i = 0), with sum_i b_i = 0. Decision variable x_ij is the flow on arc (i, j).

Minimum-cost flow — the master problem

$$ \min_{x} ; \sum_{(i,j) \in A} c_{ij}, x_{ij} $$

$$ \sum_{j:(i,j) \in A} x_{ij} ;-; \sum_{j:(j,i) \in A} x_{ji} ;=; b_i \quad \forall i \in N \qquad \text{(flow conservation)} $$

$$ 0 ;\le; x_{ij} ;\le; u_{ij} \quad \forall (i,j) \in A . $$

Every other single-commodity problem is a special case:

Problem	Specialization of min-cost flow
Shortest s-t path	b_s = 1, b_t = -1, all other b_i = 0, u = infinity
Max s-t flow	b = 0, add arc (t, s) with cost -1, u_ts = infinity, all other costs 0
Transportation	Bipartite N = sources + sinks, arcs only across
Assignment	Transportation with all supplies and demands equal to 1

Total unimodularity and free integrality

The conservation constraint matrix is the node-arc incidence matrix of a directed graph: each column has exactly one +1 and one -1. Such matrices are totally unimodular (every square submatrix has determinant in {-1, 0, +1}). By Hoffman & Kruskal (1956), a TU matrix with integral right-hand side and integral bounds defines a polyhedron with only integral vertices. Practical consequences:

• Solve the LP relaxation with simplex/barrier; any basic optimal solution is automatically an integer flow when b and u are integral. Never declare flow variables binary or integer for a pure flow problem — it only slows presolve.
• The LP duals are meaningful: the dual pi_i of node i's conservation constraint is its node potential; the reduced cost of arc (i, j) is c_ij - pi_i + pi_j. Optimality means no nonsaturated arc has negative reduced cost (Ahuja, Magnanti & Orlin 1993, Network Flows).
• Max-flow min-cut (Ford & Fulkerson 1956) is LP duality plus TU: the maximum s-t flow equals the minimum capacity over s-t cuts, and the optimal dual is the indicator of the min-cut node set.

Multicommodity flow — where integrality dies

Commodities k = 1..K, each with source s_k, sink t_k, demand d_k, and its own flow variables x^k_ij. Conservation holds per commodity; capacities couple them:

$$ \min \sum_{k} \sum_{(i,j) \in A} c_{ij}, x^{k}{ij} \qquad \text{s.t.} \qquad \sum{j} x^{k}{ij} - \sum{j} x^{k}{ji} = b^{k}{i} ;; \forall i, k, \qquad \sum_{k} x^{k}{ij} \le u{ij} ;; \forall (i,j), \qquad x \ge 0 . $$

The joint-capacity rows break total unimodularity. The splittable version is a (possibly fractional) LP. The unsplittable version — each commodity on a single path — is NP-hard already for two commodities with integral flows (Even, Itai & Shamir 1976). Expect a real integrality gap; the worked example below exhibits one (LP 11 vs unsplittable 14).

Algorithm selection

Problem	Method	Complexity	Use it when
Shortest path, c >= 0	Dijkstra, binary heap	O((n+m) log n)	Always, for nonnegative costs
Shortest path, some c < 0	Bellman-Ford-Moore (label-correcting)	O(nm)	Negative arcs or cycle detection needed
All-pairs shortest paths	Johnson (one BF reweight + n Dijkstras)	O(nm + n^2 log n)	Sparse graphs; else Floyd-Warshall O(n^3)
Max flow	Preflow-push (Goldberg & Tarjan 1988); Dinic (1970)	O(n^2 sqrt(m)); O(n^2 m)	Pure max flow / min cut
Min-cost flow	Network simplex; successive shortest path; cost scaling (Goldberg & Tarjan 1990)	Fast in practice; O(m log U) Dijkstra calls; O(nm log(nC))	Pure MCF; networkx and OR-Tools provide these
Multicommodity, splittable	LP (arc-based) or column generation (path-based)	Polynomial via LP	K and m moderate; CG when m·K explodes
Multicommodity, unsplittable	MIP, or metaheuristic over path choices	NP-hard	Exact for small K·m; heuristic beyond

Decision rule: use the dedicated algorithm while the problem stays pure; switch to the LP/MIP the moment a side constraint appears, and say explicitly in the deliverable that integrality is no longer free.

Exact Models in Gurobi

Build models with explicit constraint-builder functions: one function per constraint family, named constraints, no anonymous loops in the main routine. This makes each family unit-testable and the duals retrievable by name.

Min-cost flow LP

import gurobipy as gp
from gurobipy import GRB


def add_flow_conservation_constraints(
    model: gp.Model, x: gp.tupledict, data: dict
) -> None:
    """Add one flow-balance equality per node: outflow - inflow = balance."""
    arcs = data["arcs"]
    for i in data["nodes"]:
        out_expr = gp.quicksum(x[a] for a in arcs if a[0] == i)
        in_expr = gp.quicksum(x[a] for a in arcs if a[1] == i)
        model.addConstr(out_expr - in_expr == data["balance"][i],
                        name=f"balance[{i}]")


def add_capacity_constraints(model: gp.Model, x: gp.tupledict, data: dict) -> None:
    """Add x_a <= u_a as named constraints (instead of variable bounds) so the
    capacity dual prices are retrievable per arc."""
    for a in data["arcs"]:
        model.addConstr(x[a] <= data["capacity"][a], name=f"cap[{a[0]},{a[1]}]")


def solve_min_cost_flow(data: dict, integer: bool = False) -> dict | None:
    """Solve min-cost flow as LP (default) or IP; return flows and duals.

    integer=True is never needed for integral data (total unimodularity);
    it exists only to demonstrate that both give the same answer.
    """
    model = gp.Model("min_cost_flow")
    model.Params.OutputFlag = 0
    vtype = GRB.INTEGER if integer else GRB.CONTINUOUS
    x = model.addVars(data["arcs"], lb=0.0, vtype=vtype, name="x")
    add_flow_conservation_constraints(model, x, data)
    add_capacity_constraints(model, x, data)
    model.setObjective(
        gp.quicksum(data["cost"][a] * x[a] for a in data["arcs"]), GRB.MINIMIZE
    )
    model.optimize()
    if model.Status != GRB.OPTIMAL:
        return None
    result: dict = {
        "objective": model.ObjVal,
        "flow": {a: x[a].X for a in data["arcs"]},
    }
    if not integer:  # node potentials = duals of the conservation rows
        result["potential"] = {
            i: model.getConstrByName(f"balance[{i}]").Pi for i in data["nodes"]
        }
    return result


data = {
    "nodes": [0, 1, 2, 3],
    "arcs": [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3)],
    "balance": {0: 4, 1: 0, 2: 0, 3: -4},
    "capacity": {(0, 1): 3, (0, 2): 2, (1, 2): 2, (1, 3): 2, (2, 3): 3},
    "cost": {(0, 1): 2, (0, 2): 4, (1, 2): 1, (1, 3): 6, (2, 3): 3},
}
result = solve_min_cost_flow(data)
print(result["objective"], result["flow"])
# Expected: objective 27.0 with integral flows even though the model is an LP:
# x[0,1]=3, x[0,2]=1, x[1,2]=2, x[1,3]=1, x[2,3]=3 (total unimodularity).

Read the node potentials from result["potential"] and verify reduced-cost optimality by hand on small instances: every arc with flow strictly between 0 and u must satisfy c_ij - pi_i + pi_j = 0; every empty arc must have nonnegative reduced cost. This is the single most useful debugging check for flow LPs.

Multicommodity flow: splittable LP and unsplittable MIP

The same builder pattern, with one conservation family per (commodity, node) and one coupling family per arc. The unsplittable flag swaps continuous per-commodity flows for binary arc-usage variables scaled by demand, which forces each commodity onto a single path.

import gurobipy as gp
from gurobipy import GRB


def add_commodity_conservation_constraints(
    model: gp.Model, x: dict, data: dict
) -> None:
    """Flow balance per (commodity, node): each commodity routes its demand."""
    for k, (s_k, t_k, d_k) in enumerate(data["commodities"]):
        for i in data["nodes"]:
            rhs = d_k if i == s_k else (-d_k if i == t_k else 0.0)
            out_expr = gp.quicksum(x[k, a] for a in data["arcs"] if a[0] == i)
            in_expr = gp.quicksum(x[k, a] for a in data["arcs"] if a[1] == i)
            model.addConstr(out_expr - in_expr == rhs, name=f"bal[{k},{i}]")


def add_joint_capacity_constraints(model: gp.Model, x: dict, data: dict) -> None:
    """Couple commodities: total flow on each arc at most its capacity.
    These are the rows that destroy total unimodularity."""
    n_k = len(data["commodities"])
    for a in data["arcs"]:
        model.addConstr(
            gp.quicksum(x[k, a] for k in range(n_k)) <= data["capacity"][a],
            name=f"cap[{a[0]},{a[1]}]",
        )


def solve_multicommodity(data: dict, unsplittable: bool = False,
                         time_limit: float = 60.0) -> dict | None:
    """Arc-based multicommodity flow: LP if splittable, binary MIP if not."""
    model = gp.Model("multicommodity")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    n_k = len(data["commodities"])
    keys = [(k, a) for k in range(n_k) for a in data["arcs"]]
    if unsplittable:
        y = model.addVars(keys, vtype=GRB.BINARY, name="y")
        # x[k,a] = d_k * y[k,a]: commodity k uses arc a fully or not at all.
        x = {(k, a): data["commodities"][k][2] * y[k, a] for (k, a) in keys}
    else:
        x = model.addVars(keys, lb=0.0, name="x")
    add_commodity_conservation_constraints(model, x, data)
    add_joint_capacity_constraints(model, x, data)
    model.setObjective(
        gp.quicksum(data["cost"][a] * x[k, a] for (k, a) in keys), GRB.MINIMIZE
    )
    model.optimize()
    ok = model.Status == GRB.OPTIMAL or (
        model.Status == GRB.TIME_LIMIT and model.SolCount > 0
    )
    if not ok:
        return None
    if unsplittable:
        flow = {ka: data["commodities"][ka[0]][2] * y[ka].X for ka in keys}
    else:
        flow = {ka: x[ka].X for ka in keys}
    return {"objective": model.ObjVal, "flow": flow, "gap": model.MIPGap
            if unsplittable else 0.0}


data = {
    "nodes": [0, 1, 2, 3],
    "arcs": [(0, 2), (1, 2), (2, 3), (0, 3), (1, 3)],
    "capacity": {(0, 2): 3, (1, 2): 3, (2, 3): 3, (0, 3): 3, (1, 3): 3},
    "cost": {(0, 2): 1, (1, 2): 1, (2, 3): 1, (0, 3): 5, (1, 3): 5},
    "commodities": [(0, 3, 2.0), (1, 3, 2.0)],  # (source, sink, demand)
}
lp = solve_multicommodity(data)
ip = solve_multicommodity(data, unsplittable=True)
print(lp["objective"], ip["objective"])
# Expected: splittable LP objective 11.0 (3 units share the cheap hub arc
# (2,3), 1 unit takes a direct arc); unsplittable objective 14.0 (one
# commodity hub, one direct). The gap exists because the joint-capacity
# matrix is not totally unimodular.

For large instances the arc-based LP has m·K variables and grows fast; switch to the path-based formulation with column generation (see Advanced Techniques) before fighting the LP size.

Dedicated Algorithms: Shortest Paths and networkx

Use these when the problem is pure. They are 10-1000x faster than the LP route and need no license.

Dijkstra — label-setting, nonnegative costs only

Dijkstra (1959) permanently settles nodes in order of distance. Correctness relies on c >= 0: once a node leaves the heap with its final label, no later path can improve it. With a binary heap and lazy deletion it runs in O((n + m) log n).

import heapq


def dijkstra(
    adj: dict[int, list[tuple[int, float]]], source: int
) -> tuple[dict[int, float], dict[int, int]]:
    """Label-setting shortest paths from source; requires c >= 0 on all arcs.

    adj[i] = list of (j, cost_ij). Returns (distance, predecessor) maps.
    Lazy deletion: stale heap entries are skipped via the settled set.
    """
    dist: dict[int, float] = {source: 0.0}
    pred: dict[int, int] = {}
    heap: list[tuple[float, int]] = [(0.0, source)]
    settled: set[int] = set()
    while heap:
        d_i, i = heapq.heappop(heap)
        if i in settled:
            continue
        settled.add(i)
        for j, c_ij in adj.get(i, []):
            d_j = d_i + c_ij
            if d_j < dist.get(j, float("inf")):
                dist[j] = d_j
                pred[j] = i
                heapq.heappush(heap, (d_j, j))
    return dist, pred


adj = {0: [(1, 2.0), (2, 4.0)], 1: [(2, 1.0), (3, 6.0)], 2: [(3, 3.0)]}
dist, pred = dijkstra(adj, 0)
print(dist[3], pred)
# Expected: dist[3] == 6.0 along 0 -> 1 -> 2 -> 3; pred == {1: 0, 2: 1, 3: 2}.

Bellman-Ford-Moore — label-correcting, handles negative costs

The queue-based form below is the classic label-correcting method (Bellman 1958; Ford 1956; Moore 1959): node labels stay tentative and may improve many times. It tolerates negative arc costs and proves negative cycles. Variants differ only in queue discipline — FIFO here; deque insertion (d'Esopo-Pape, Pape 1974) and small-label-first are drop-in replacements that often run faster on road-like networks but have bad worst cases.

from collections import deque


def bellman_ford(
    arcs: list[tuple[int, int, float]], n_nodes: int, source: int
) -> tuple[list[float], list[int], bool]:
    """FIFO label-correcting shortest paths; allows negative arc costs.

    Returns (dist, pred, has_negative_cycle). A node dequeued more than
    n_nodes times proves a negative cycle reachable from the source.
    O(n*m) worst case; close to O(m) on sparse practical graphs.
    """
    inf = float("inf")
    dist = [inf] * n_nodes
    pred = [-1] * n_nodes
    dist[source] = 0.0
    out: dict[int, list[tuple[int, float]]] = {}
    for i, j, c in arcs:
        out.setdefault(i, []).append((j, c))
    in_queue = [False] * n_nodes
    dequeues = [0] * n_nodes
    queue: deque[int] = deque([source])
    in_queue[source] = True
    while queue:
        i = queue.popleft()
        in_queue[i] = False
        dequeues[i] += 1
        if dequeues[i] > n_nodes:
            return dist, pred, True
        for j, c in out.get(i, []):
            if dist[i] + c < dist[j] - 1e-12:
                dist[j] = dist[i] + c
                pred[j] = i
                if not in_queue[j]:
                    queue.append(j)
                    in_queue[j] = True
    return dist, pred, False


arcs = [(0, 1, 4.0), (0, 2, 2.0), (2, 1, -3.0), (1, 3, 2.0), (2, 3, 4.0)]
dist, pred, neg = bellman_ford(arcs, 4, 0)
print(dist[3], neg)
# Expected: dist[3] == 1.0 via 0 -> 2 -> 1 -> 3 (2 - 3 + 2); neg == False.

For repeated shortest-path calls on a graph with some negative costs, run Bellman-Ford once, reweight every arc to c_ij - pi_i + pi_j >= 0 with the resulting potentials, then use Dijkstra for all later queries (Johnson 1977). This is exactly how successive-shortest-path min-cost flow stays fast.

networkx for max flow, min cut, and min-cost flow

import networkx as nx


def networkx_flows() -> None:
    """Reference solves: max-flow/min-cut and min-cost flow with networkx."""
    # Max flow / min cut (preflow-push is the default algorithm).
    g = nx.DiGraph()
    for (i, j), u in {(0, 1): 3, (0, 2): 2, (1, 2): 1,
                      (1, 3): 2, (2, 3): 3}.items():
        g.add_edge(i, j, capacity=u)
    value, flow = nx.maximum_flow(g, 0, 3)
    cut_value, (side_s, side_t) = nx.minimum_cut(g, 0, 3)
    print(value, cut_value, sorted(side_s))
    # Expected: value == 5, cut_value == 5, side_s == [0]
    # (the min cut is {(0,1),(0,2)} with capacity 3 + 2 = 5).

    # Min-cost flow via network simplex. Convention: node attribute
    # 'demand' is negative for supply. Data must be integral.
    h = nx.DiGraph()
    h.add_node(0, demand=-4)
    h.add_node(3, demand=4)
    arc_data = {(0, 1): (3, 2), (0, 2): (2, 4), (1, 2): (2, 1),
                (1, 3): (2, 6), (2, 3): (3, 3)}
    for (i, j), (u, c) in arc_data.items():
        h.add_edge(i, j, capacity=u, weight=c)
    flow_dict = nx.min_cost_flow(h)
    print(nx.cost_of_flow(h, flow_dict))
    # Expected: 27 — identical to the Gurobi LP on the same instance.


networkx_flows()

networkx is the right prototyping tool up to roughly 10^5 arcs. Beyond that, use OR-Tools (ortools.graph.python.min_cost_flow and max_flow), which wrap C++ cost-scaling and push-relabel implementations (Goldberg 1997, "An efficient implementation of a scaling minimum-cost flow algorithm").

Metaheuristic: GA for Unsplittable Multicommodity Flow

Pure single-commodity flow problems never need a metaheuristic — the polynomial algorithms above are exact and fast. The NP-hard member of this family is unsplittable multicommodity flow, and a genetic algorithm over path-index chromosomes matches it well: precompute a small candidate-path set per commodity (k-shortest paths via networkx.shortest_simple_paths, or repeated Dijkstra with perturbed arc costs), then search over the index vector that picks one path per commodity. Capacity violations are handled by a static penalty. Operator and parameter details are covered in the genetic-algorithms and constraint-handling-techniques skills; the encoding here keeps everything vectorizable.

import numpy as np


def path_arc_matrices(
    candidates: list[list[list[tuple[int, int]]]],
    arcs: list[tuple[int, int]],
    demands: np.ndarray,
    costs: dict[tuple[int, int], float],
) -> tuple[np.ndarray, np.ndarray]:
    """Build load tensor L[k,p,a] (= d_k if commodity k's path p uses arc a)
    and path-cost matrix C[k,p]; unused path slots get cost +inf."""
    n_k, n_p = len(candidates), max(len(c) for c in candidates)
    arc_idx = {a: t for t, a in enumerate(arcs)}
    load = np.zeros((n_k, n_p, len(arcs)))
    cost = np.full((n_k, n_p), np.inf)
    for k, paths in enumerate(candidates):
        for p, path in enumerate(paths):
            load[k, p, [arc_idx[a] for a in path]] = demands[k]
            cost[k, p] = demands[k] * sum(costs[a] for a in path)
    return load, cost


def ga_unsplittable_mcf(
    load: np.ndarray, path_cost: np.ndarray, capacity: np.ndarray,
    n_paths: np.ndarray, penalty: float = 1_000.0, pop: int = 60,
    gens: int = 200, seed: int = 0,
) -> tuple[np.ndarray, float]:
    """Vectorized GA over path-index chromosomes (one gene per commodity).
    See the genetic-algorithms skill for operator design rationale."""
    rng = np.random.default_rng(seed)
    n_k = load.shape[0]
    k_idx = np.arange(n_k)

    def fitness(pop_mat: np.ndarray) -> np.ndarray:
        arc_load = load[k_idx, pop_mat, :].sum(axis=1)        # (pop, arcs)
        excess = np.maximum(arc_load - capacity, 0.0).sum(axis=1)
        return path_cost[k_idx, pop_mat].sum(axis=1) + penalty * excess

    P = rng.integers(0, n_paths, size=(pop, n_k))
    fit = fitness(P)
    for _ in range(gens):
        a = rng.integers(0, pop, (2, pop))                    # tournaments
        b = rng.integers(0, pop, (2, pop))
        pa = np.where(fit[a[0]] < fit[a[1]], a[0], a[1])
        pb = np.where(fit[b[0]] < fit[b[1]], b[0], b[1])
        mask = rng.random((pop, n_k)) < 0.5                   # uniform xover
        children = np.where(mask, P[pa], P[pb])
        mut = rng.random((pop, n_k)) < 1.0 / n_k              # random reset
        children = np.where(mut, rng.integers(0, n_paths, (pop, n_k)), children)
        merged = np.vstack([P, children])                     # (mu + lambda)
        merged_fit = np.concatenate([fit, fitness(children)])
        keep = np.argpartition(merged_fit, pop - 1)[:pop]
        P, fit = merged[keep], merged_fit[keep]
    best = int(np.argmin(fit))
    return P[best], float(fit[best])


arcs = [(0, 2), (1, 2), (2, 3), (0, 3), (1, 3)]
costs = {(0, 2): 1.0, (1, 2): 1.0, (2, 3): 1.0, (0, 3): 5.0, (1, 3): 5.0}
capacity = np.array([3.0, 3.0, 3.0, 3.0, 3.0])
demands = np.array([2.0, 2.0])
candidates = [
    [[(0, 2), (2, 3)], [(0, 3)]],   # commodity 0: hub path or direct arc
    [[(1, 2), (2, 3)], [(1, 3)]],   # commodity 1: hub path or direct arc
]
load, path_cost = path_arc_matrices(candidates, arcs, demands, costs)
chrom, best_fit = ga_unsplittable_mcf(load, path_cost, capacity,
                                      np.array([2, 2]), seed=1)
print(chrom, best_fit)
# Expected: best_fit == 14.0 — one commodity through the hub, the other
# direct; equals the unsplittable Gurobi optimum on the same instance.

Two design notes. First, the candidate-path set bounds solution quality: if the optimal path of some commodity is missing, the GA cannot find it, so generate 5-20 diverse candidates per commodity and report which were used. Second, always benchmark the GA against the splittable LP value — it is a valid lower bound, and (GA - LP) / LP is an honest quality gap even when the MIP is too large to solve.

Instance Generation and Validation

Seeded generators (feasible by construction)

Random balances with random capacities are usually infeasible. Generate instances by routing a known flow first and sizing capacities to cover it: feasibility is then guaranteed, and the routed flow gives a built-in upper bound on the optimal cost.

import numpy as np


def generate_mcf_instance(
    n_nodes: int, arc_density: float, n_supply: int, seed: int
) -> dict:
    """Random feasible min-cost flow instance.

    Builds a DAG skeleton (arcs i -> j only for i < j) including the chain
    0 -> 1 -> ... -> n-1 for connectivity, routes random supplies to the last
    node along node-increasing walks, then sets capacities to cover the
    routed flow plus slack. Feasible by construction.
    """
    rng = np.random.default_rng(seed)
    nodes = list(range(n_nodes))
    arc_set = {(i, i + 1) for i in range(n_nodes - 1)}
    for i in range(n_nodes):
        for j in range(i + 2, n_nodes):
            if rng.random() < arc_density:
                arc_set.add((i, j))
    arcs = sorted(arc_set)
    out: dict[int, list[int]] = {}
    for i, j in arcs:
        out.setdefault(i, []).append(j)
    used = {a: 0.0 for a in arcs}
    balance = {i: 0.0 for i in nodes}
    sink = n_nodes - 1
    for _ in range(n_supply):
        amount = float(rng.integers(1, 10))
        i = int(rng.integers(0, sink))
        balance[i] += amount
        while i != sink:                      # node-increasing walk to sink
            j = int(rng.choice(out[i]))
            used[(i, j)] += amount
            i = j
        balance[sink] -= amount
    capacity = {a: used[a] + float(rng.integers(0, 8)) for a in arcs}
    cost = {a: float(rng.integers(1, 20)) for a in arcs}
    return {"nodes": nodes, "arcs": arcs, "balance": balance,
            "capacity": capacity, "cost": cost}


def generate_multicommodity_instance(
    n_nodes: int, arc_density: float, n_commodities: int, seed: int
) -> dict:
    """Feasible multicommodity instance on the same DAG skeleton: each
    commodity gets a routed witness path; joint capacities cover the sum."""
    rng = np.random.default_rng(seed)
    base = generate_mcf_instance(n_nodes, arc_density, n_supply=1, seed=seed)
    arcs, out = base["arcs"], {}
    for i, j in arcs:
        out.setdefault(i, []).append(j)
    used = {a: 0.0 for a in arcs}
    commodities: list[tuple[int, int, float]] = []
    for _ in range(n_commodities):
        s = int(rng.integers(0, n_nodes - 1))
        d = float(rng.integers(1, 6))
        i = s
        while i != n_nodes - 1:
            j = int(rng.choice(out[i]))
            used[(i, j)] += d
            i = j
        commodities.append((s, n_nodes - 1, d))
    capacity = {a: used[a] + float(rng.integers(0, 5)) for a in arcs}
    return {"nodes": base["nodes"], "arcs": arcs, "capacity": capacity,
            "cost": base["cost"], "commodities": commodities}


inst = generate_mcf_instance(n_nodes=8, arc_density=0.3, n_supply=5, seed=42)
print(len(inst["arcs"]), sum(inst["balance"].values()))
# Expected: a connected DAG with the printed arc count for seed 42, and
# balances summing to exactly 0.0 (feasible by construction).

Report the generator parameters (n_nodes, density, seed) with every experiment so instances are reproducible; vary seeds, never instances by hand.

Independent validator

Never trust the model that produced the solution. The validator recomputes feasibility and objective from raw instance data only.

def validate_flow(
    data: dict, flow: dict, tol: float = 1e-6
) -> tuple[bool, float, list[str]]:
    """Independent check of a single-commodity flow: sign, capacity,
    conservation at every node, and recomputed objective."""
    errors: list[str] = []
    for a in data["arcs"]:
        f = flow.get(a, 0.0)
        if f < -tol:
            errors.append(f"negative flow {f:.6g} on arc {a}")
        if f > data["capacity"][a] + tol:
            errors.append(
                f"capacity exceeded on {a}: {f:.6g} > {data['capacity'][a]}")
    for i in data["nodes"]:
        net = (sum(flow.get(a, 0.0) for a in data["arcs"] if a[0] == i)
               - sum(flow.get(a, 0.0) for a in data["arcs"] if a[1] == i))
        if abs(net - data["balance"][i]) > tol:
            errors.append(
                f"conservation violated at node {i}: "
                f"net {net:.6g} != balance {data['balance'][i]}")
    objective = sum(data["cost"][a] * flow.get(a, 0.0) for a in data["arcs"])
    return len(errors) == 0, objective, errors


def validate_multicommodity(
    data: dict, flow: dict, tol: float = 1e-6
) -> tuple[bool, float, list[str]]:
    """Check per-commodity conservation plus joint capacities; recompute cost.
    flow is keyed by (k, arc)."""
    errors: list[str] = []
    n_k = len(data["commodities"])
    for k, (s_k, t_k, d_k) in enumerate(data["commodities"]):
        for i in data["nodes"]:
            rhs = d_k if i == s_k else (-d_k if i == t_k else 0.0)
            net = (sum(flow.get((k, a), 0.0)
                       for a in data["arcs"] if a[0] == i)
                   - sum(flow.get((k, a), 0.0)
                         for a in data["arcs"] if a[1] == i))
            if abs(net - rhs) > tol:
                errors.append(f"commodity {k} conservation broken at node {i}")
    for a in data["arcs"]:
        total = sum(flow.get((k, a), 0.0) for k in range(n_k))
        if total > data["capacity"][a] + tol:
            errors.append(f"joint capacity exceeded on arc {a}: {total:.6g}")
    objective = sum(data["cost"][a] * flow.get((k, a), 0.0)
                    for k in range(n_k) for a in data["arcs"])
    return len(errors) == 0, objective, errors


data = {
    "nodes": [0, 1, 2, 3],
    "arcs": [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3)],
    "balance": {0: 4, 1: 0, 2: 0, 3: -4},
    "capacity": {(0, 1): 3, (0, 2): 2, (1, 2): 2, (1, 3): 2, (2, 3): 3},
    "cost": {(0, 1): 2, (0, 2): 4, (1, 2): 1, (1, 3): 6, (2, 3): 3},
}
flow = {(0, 1): 3.0, (0, 2): 1.0, (1, 2): 2.0, (1, 3): 1.0, (2, 3): 3.0}
print(validate_flow(data, flow))
# Expected: (True, 27.0, []) — confirms the Gurobi and networkx results.

Wire the validator into every solve path: exact, heuristic, and any solution read back from a file. Cross-check exact vs heuristic on small instances — on the multicommodity example above, the validator must report 14.0 for both the MIP solution and the GA solution.

Advanced Techniques

Node potentials, reduced costs, and warm-started pricing

The dual pi_i of node i turns arc costs into reduced costs c^pi_ij = c_ij - pi_i + pi_j. Successive-shortest-path min-cost flow keeps all reduced costs nonnegative, so each augmentation step is a Dijkstra call even when original costs are negative (Ahuja, Magnanti & Orlin 1993, ch. 9). When flow problems sit inside a loop — column-generation pricing, Benders subproblems, rolling-horizon resolves — keep the potentials from the previous solve as the starting reweighting: the first Dijkstra after a small data change then touches only the affected region of the graph.

Scaling algorithms for large instances

Successive shortest path alone is pseudo-polynomial (one augmentation per unit of flow in the worst case). Capacity scaling fixes this by augmenting only along paths with residual capacity at least Delta, halving Delta over O(log U) phases, giving O(m log U) Dijkstra calls. Cost scaling (Goldberg & Tarjan 1990) refines an epsilon-optimality condition on reduced costs and runs in O(n m log(nC)); it is the algorithm behind OR-Tools' min-cost flow. Network simplex — simplex specialized to spanning-tree bases with O(1) pivot pricing per arc — is exponential in theory but typically the fastest exact method in practice, and Orlin (1997) gives a polynomial variant. Practical order of attack: network simplex (networkx, Gurobi Method=0 on the LP) first, cost scaling (OR-Tools) when m exceeds ~10^5.

Lagrangian relaxation of joint capacities

For multicommodity flow, dualize the joint-capacity rows with multipliers w_ij >= 0. The Lagrangian subproblem separates into K independent shortest paths under costs c_ij + w_ij — each solved by Dijkstra in milliseconds. Update w by subgradient steps on the capacity violations. Because the subproblems have the integrality property, the best Lagrangian bound equals the LP relaxation value (Geoffrion 1974, "Lagrangean relaxation for integer programming") — you gain speed and a primal heuristic (route commodities greedily under the final w), not a tighter bound.

Path-based formulation and column generation

Replace arc variables with path variables f_p >= 0 for p in P_k, the set of s_k-t_k paths. The master LP has one convexity row per commodity and one capacity row per arc; the pricing problem for commodity k is a shortest path under reduced costs c_ij - w_ij (w = capacity duals), so new columns come from Dijkstra. This scales far beyond the arc formulation when K is large and is the standard route to exact unsplittable solutions via branch-and-price (Barnhart, Hane & Vance 2000, "Using branch-and-price-and- cut to solve origin-destination integer multicommodity flow problems").

Graph transformations catalog

Reduce nonstandard features to the standard model instead of writing new algorithms:

Feature	Transformation
Node capacity / node cost at i	Split i into i_in, i_out joined by one arc carrying the capacity/cost
Arc lower bound l_ij > 0	Send l_ij permanently: set u' = u - l, add l_ij to b_j, subtract from b_i
Undirected edge, c >= 0	Two antiparallel arcs with the same capacity and cost
Multiple sources/sinks (max flow)	Super-source/super-sink with infinite-capacity arcs
Max flow as min-cost flow	Zero costs, add arc (t, s) with cost -1 and infinite capacity
Unbalanced supply/demand	Dummy node absorbing the surplus via zero-cost arcs

Practical Challenges

Supplies and demands do not balance. nx.min_cost_flow raises NetworkXUnfeasible and the Gurobi LP reports infeasible when sum b_i != 0. Decide what surplus means in the application — lost sales, inventory, overtime — and add a dummy node that absorbs it through arcs whose cost prices that interpretation. A zero-cost dummy silently hides modeling errors; give it a visible cost.

networkx min-cost flow with float data. nx.min_cost_flow and nx.network_simplex assume integer capacities, demands, and weights; float input can produce wrong answers or non-termination. Scale to integers (multiply by 10^k, round, record k) or use the Gurobi/HiGHS LP, which is float-native.

Dijkstra silently wrong with negative arcs. Dijkstra never revisits a settled node, so one negative arc can make it return a non-shortest path with no error raised. Scan min(c) at the boundary of your code; dispatch to Bellman-Ford or Johnson reweighting when it is negative. Add a regression test with a known negative-arc instance.

Infeasible instance with no explanation. Conservation infeasibility is a cut problem: build the super-source/super-sink max-flow instance; if the max flow is less than total supply, the min cut returned by nx.minimum_cut names exactly the bottleneck arc set separating supply from demand. Report that cut, not just "infeasible". On the LP route, model.computeIIS() gives the conflicting constraint set.

Parallel arcs disappear. nx.DiGraph keeps one edge per node pair, so two cables between the same substations merge silently. Use nx.MultiDiGraph (flow functions accept it) or key arcs as (i, j, idx) in your own models and generators.

Side constraints quietly destroy integrality. Adding a budget row, fixed charges, or arc-choice binaries to a flow LP breaks total unimodularity; solutions may turn fractional and rounding them violates conservation. The moment a non-network row enters, declare the integer variables integer, treat the model as a true MIP, and re-estimate solve time — see milp-modeling-gurobi.

Python loops too slow on large graphs. Dict-based Dijkstra at 10^6 arcs costs seconds per call; inside a pricing loop that is fatal. Move to scipy.sparse.csgraph (CSR adjacency, C speed) or OR-Tools. Keep the pure- Python version as the reference implementation for tests.

Degenerate duals confuse sensitivity stories. Flow LPs are massively degenerate; potentials are not unique, and two solvers can return different pi vectors for the same optimum. Only reduced-cost signs and objective ranges are reliable claims; never report a single dual vector as "the" marginal prices without checking uniqueness on a perturbed instance.

Tools & Libraries

Library	When to use	Note
networkx	Prototyping, graphs up to ~10^5 arcs	maximum_flow, minimum_cut, min_cost_flow, network_simplex; integer data required for MCF
gurobipy	Flow + side constraints, multicommodity, duals needed	LP/MIP route; TU gives integral LP solutions for pure flow
OR-Tools graph	Production scale, 10^6+ arcs, no license	min_cost_flow (cost scaling), max_flow (push-relabel); C++ speed, integer data
scipy.sparse.csgraph	Many shortest-path calls on a fixed graph	dijkstra, bellman_ford, johnson, plus maximum_flow in scipy.sparse
HiGHS (via scipy.optimize.linprog or highspy)	LP route without a Gurobi license	Handles the arc-based formulations in this skill
igraph / graph-tool	Very large graphs, analysis-heavy workflows	C-backed; fewer flow-specific APIs than OR-Tools

Output Format

A complete network-flow deliverable contains:

1. Problem classification. One line: which flow problem, single or multicommodity, pure or with side constraints, and the chosen method with the reason ("pure MCF, integral data -> network simplex via networkx").

2. Model summary table.

   Item	Value
   Nodes / arcs / commodities	n, m, K
   Variables / constraints (LP or MIP route)	counts
   Integrality	free (TU) or enforced (MIP)
   Solver and algorithm	e.g. Gurobi LP dual simplex; OR-Tools cost scaling

3. Solution report. Objective value; table of nonzero arc flows (arc, flow, capacity, cost, reduced cost); node potentials when the LP route was used; the min-cut arc set when max flow was solved; for MIPs, the bound and final gap; runtime and seed.

4. Validation line. Output of the independent validator on the reported flow: feasible=True, objective=..., 0 violations. A deliverable without this line is incomplete.

5. Heuristic extras (when a metaheuristic was used). Candidate-set size per commodity, generations, population, penalty value, best/mean over seeds, and the gap to the splittable-LP lower bound.

6. Artifacts. solution.csv (columns: commodity, tail, head, flow), instance.json (generator name + parameters + seed), and the run log.

Questions to Ask

• Is this a pure flow problem, or are there side constraints (fixed charges, budgets, logical rules) that break the network structure?
• One commodity or several? If several: may a commodity split across paths?
• Are capacities, demands, and costs integral, or floats that need scaling?
• Any negative arc costs — and could there be a negative cycle?
• Directed arcs, undirected edges, or a mix?
• How large: nodes, arcs, commodities? One solve or thousands in a loop?
• Do supplies and demands balance, and what should happen to any surplus?
• Are duals, sensitivity ranges, or a min-cut certificate part of the deliverable, or only the flow itself?
• Is a Gurobi license available, or should everything run on open tools?
• What time budget separates "exact MIP" from "heuristic with LP bound"?

Related Skills

• linear-programming-fundamentals — when you need the duality, reduced- cost, and sensitivity background behind node potentials and total unimodularity.
• assignment-problems — when the network is bipartite with unit supplies and demands; assignment is the min-cost flow special case with dedicated algorithms.
• milp-modeling-gurobi — when side constraints break the network structure and a general MIP must be built and solved properly.
• dynamic-programming — when shortest paths carry resource constraints; labeling algorithms are DP over the node-resource state space and feed pricing loops.