Assignment Problems

Assignment Problems

You are an expert in assignment problems: linear assignment (LAP), generalized assignment (GAP), and bottleneck assignment. This skill covers exact methods (Hungarian algorithm, LP with total unimodularity, MIP), Lagrangian bounds, and local-search heuristics, plus instance generation and independent solution validation. Use the framework below to classify the variant, pick the cheapest adequate method, implement it, and verify the result.

Initial Assessment

Establish these points before formulating or recommending a method:

• Cardinality structure. One agent per task and one task per agent

(one-to-one, LAP)? Or can one agent take several tasks subject to a capacity (many-to-one, GAP)? This single distinction separates a polynomial problem from an NP-hard one.

• Objective sense and shape. Minimize total cost, maximize total profit,

or minimize the worst single cost (bottleneck)? Mixed conventions are the most common source of wrong answers; fix the sense first.

• Sizes. Number of agents m, tasks n. LAP with n up to ~10,000 is

routine for scipy.optimize.linear_sum_assignment. GAP with m·n up to ~10^5 binaries is usually fine for a MIP solver; beyond that plan for Lagrangian bounds plus a heuristic.

• Balanced or rectangular. Equal numbers on both sides? If not, decide

whether unmatched rows/columns are allowed and what they cost.

• Forbidden pairs. Are some (agent, task) combinations disallowed? Plan

to encode them as np.inf (scipy) or by omitting variables (MIP), not as fragile big-M costs.

• Resource data type (GAP). Integer resource consumptions enable

knapsack DP in the Lagrangian subproblems; float data must be scaled or the subproblems solved as small MIPs.

• Solve count. One-off solve, or LAP/GAP called thousands of times

inside a heuristic or branch-and-bound loop? The embedded case changes the tooling (warm starts, C-backed LAP libraries, candidate lists).

• Solver availability. Gurobi license present? If not, scipy +

open-source MIP (HiGHS) covers everything in this skill.

• Optimality requirement. LAP and bottleneck are always exact. For GAP,

ask whether a proven optimum is required or a bounded-gap heuristic solution within a time budget suffices.

• Validation path. Agree up front that every reported solution passes an

independent feasibility and objective check (provided below).

Problem Variants and Formulation

Linear assignment problem (LAP)

Given an n×n cost matrix C = (c_ij), choose a permutation assigning each row to exactly one column:

$$ \min \sum{i=1}^{n}\sum{j=1}^{n} c{ij}\, x{ij} \quad \text{s.t.} \quad \sum{j} x{ij} = 1 \;\; \forall i, \qquad \sum{i} x{ij} = 1 \;\; \forall j, \qquad x_{ij} \ge 0 . $$

No integrality constraints are written, and none are needed. The constraint matrix is the node-edge incidence matrix of a bipartite graph, which is totally unimodular: every square submatrix has determinant in {-1, 0, +1}, so every basic feasible solution of the LP is integral. Equivalently, the extreme points of the feasible polytope (the Birkhoff polytope of doubly stochastic matrices) are exactly the permutation matrices (Birkhoff 1946). The Hungarian method (Kuhn 1955; Munkres 1957) is therefore a primal-dual LP algorithm: it keeps dual potentials ui, vj with ui + vj ≤ c_ij and grows a matching on tight edges; at optimality Σu + Σv equals the optimal cost. Modern dense implementations follow Jonker & Volgenant (1987), O(n³).

Bottleneck assignment problem

Same feasible set, different objective: minimize the largest cost used,

$$ \min{x} \; \max \{\, c{ij} : x_{ij} = 1 \,\}. $$

Use it when the slowest pairing determines system performance (parallel workers, latest completion time, worst-case latency). Threshold search solves it exactly: the optimum is one of the O(n²) distinct cost values, and a threshold t is feasible iff the edges with c_ij ≤ t admit a perfect matching. Hopcroft-Karp gives ~O(n^2.5 log n); faster methods exist (Gabow & Tarjan 1988).

Generalized assignment problem (GAP)

Agents I = {1,…,m} with capacities bi; jobs J = {1,…,n}. Agent i spends resource aij and cost c_ij on job j:

$$ \min \sum{i \in I}\sum{j \in J} c{ij}\, x{ij} $$

$$ \sum{i \in I} x{ij} = 1 \quad \forall j \in J \qquad \text{(each job done by exactly one agent)} $$

$$ \sum{j \in J} a{ij}\, x{ij} \le bi \quad \forall i \in I \qquad \text{(agent capacity)}, \qquad x_{ij} \in \{0,1\}. $$

The capacity constraints destroy total unimodularity: the LP relaxation is fractional, the problem is NP-hard, and even deciding feasibility is NP-complete (Martello & Toth 1990, Knapsack Problems, ch. 7). Much of the literature states GAP as profit maximization; convert with c'ij = maxkl(ckl) − cij and keep one convention throughout your code. LAP is the special case m = n, aij = 1, bi = 1 with equality capacities.

Method selection

| Situation | Method | |---|---| | One-to-one, dense costs, n ≤ ~10^4 | scipy.optimize.linear_sum_assignment | | One-to-one, sparse costs | scipy.sparse.csgraph.min_weight_full_bipartite_matching | | One-to-one, need duals / sensitivity | LP in gurobipy, read Pi (or own Hungarian potentials) | | Min-max fairness objective | Bottleneck threshold search | | Many-to-one with capacities, moderate size | GAP MIP in gurobipy | | GAP, large or time-boxed | Lagrangian bound + local search / ejection chains | | Costs depend on pairs of assignments | Not LAP/GAP — see quadratic-assignment-problem | | Assignment is a subproblem of a flow network | See network-flow-optimization (min-cost flow) |

Linear and Bottleneck Assignment

scipy: the default LAP tool

linear_sum_assignment implements a Jonker-Volgenant-type shortest augmenting path algorithm in C. It handles rectangular matrices (every row matched when rows ≤ columns), maximization, and np.inf for forbidden pairs.

import numpy as np
from scipy.optimize import linear_sum_assignment

def solve_lap(cost: np.ndarray, maximize: bool = False) -> tuple[np.ndarray, np.ndarray, float]:
    """Solve a (possibly rectangular) LAP with scipy's JV-type solver.

    Returns (rows, cols, total); row rows[k] is matched to column cols[k].
    np.inf entries mark forbidden pairs; scipy raises ValueError when no
    feasible complete assignment exists.
    """
    rows, cols = linear_sum_assignment(cost, maximize=maximize)
    return rows, cols, float(cost[rows, cols].sum())

cost = np.array([[4.0, 1.0, 3.0],
                 [2.0, 0.0, 5.0],
                 [3.0, 2.0, 2.0]])
print(solve_lap(cost))
# Expected: rows [0 1 2], cols [1 0 2], total 5.0

forbidden = cost.copy()
forbidden[0, 1] = np.inf                  # row 0 may not take column 1
print(solve_lap(forbidden))
# Expected: total 6.0 (two optimal supports exist for this matrix)

tasks = np.array([[9.0, 4.0, 6.0, 2.0],
                  [3.0, 8.0, 5.0, 7.0]])  # 2 workers, 4 tasks: rectangular
rows, cols, total = solve_lap(tasks)
print(rows, cols, total)
# Expected: both rows matched to distinct columns, total 5.0 (2.0 + 3.0)

A reference Hungarian implementation

Owning a transparent O(n³) implementation is useful when you must instrument the algorithm (extract potentials, warm-start, embed where C extensions are unavailable). This is the shortest-augmenting-path variant; costs must be finite.

import numpy as np

def hungarian(cost: np.ndarray) -> tuple[np.ndarray, float]:
    """Solve a square LAP by shortest augmenting paths, O(n^3).

    Maintains dual potentials u (rows) and v (columns) with
    u[i] + v[j]  tuple[np.ndarray, float]:
    """Solve the assignment problem as a pure LP.

    Total unimodularity of the bipartite incidence matrix guarantees an
    integral optimal vertex, so binary variables are unnecessary.
    Returns (col_of_row, optimal value).
    """
    n = cost.shape[0]
    model = gp.Model("lap_lp")
    model.Params.OutputFlag = 0
    model.Params.Method = 0                  # primal simplex -> vertex solution
    x = model.addVars(n, n, lb=0.0, ub=1.0, name="x")
    model.setObjective(
        gp.quicksum(cost[i, j] * x[i, j] for i in range(n) for j in range(n)),
        GRB.MINIMIZE,
    )
    row_c = model.addConstrs((x.sum(i, "*") == 1 for i in range(n)), name="row")
    col_c = model.addConstrs((x.sum("*", j) == 1 for j in range(n)), name="col")
    model.optimize()
    if model.Status != GRB.OPTIMAL:
        raise RuntimeError(f"unexpected LP status {model.Status}")
    sol = np.array([[x[i, j].X for j in range(n)] for i in range(n)])
    assert np.allclose(sol, sol.round()), "TU guarantees an integral vertex"
    col_of_row = sol.argmax(axis=1)
    # The LP duals are Hungarian potentials: u[i] + v[j]  tuple[np.ndarray, float]:
    """Solve the bottleneck (min-max) assignment problem exactly.

    Binary search over the sorted distinct cost values; feasibility of a
    threshold t is a perfect-matching test on edges with cost = 0).all():               # perfect matching exists
            best_match, best_val = match.copy(), float(values[mid])
            hi = mid - 1
        else:
            lo = mid + 1
    if best_match is None:
        raise ValueError("no perfect matching exists (forbidden pairs block all)")
    return best_match.astype(int), best_val

cost = np.array([[3.0, 8.0, 6.0],
                 [2.0, 4.0, 9.0],
                 [7.0, 5.0, 1.0]])
col_of_row, value = bottleneck_assignment(cost)
print(col_of_row, value)
# Expected: col_of_row [0 1 2], bottleneck value 4.0
# (threshold 3 fails: rows 0 and 1 would both need column 0)

Generalized Assignment: Exact MIP, Instances, Validation

Exact MIP with explicit constraint builders

Keep each constraint family in its own named builder function. This makes the model auditable, testable in isolation, and easy to extend (e.g., adding assignment restrictions or multiple resources later).

import gurobipy as gp
from gurobipy import GRB
import numpy as np

def add_job_assignment_constraints(model: gp.Model, x: gp.tupledict,
                                   n_agents: int, n_jobs: int) -> None:
    """Every job is assigned to exactly one agent."""
    model.addConstrs(
        (gp.quicksum(x[i, j] for i in range(n_agents)) == 1 for j in range(n_jobs)),
        name="assign",
    )

def add_capacity_constraints(model: gp.Model, x: gp.tupledict,
                             a: np.ndarray, b: np.ndarray) -> None:
    """Resource consumed on each agent stays within its capacity."""
    m, n = a.shape
    model.addConstrs(
        (gp.quicksum(a[i, j] * x[i, j] for j in range(n))  tuple[gp.Model, gp.tupledict]:
    """Assemble the min-cost GAP MIP from its constraint builders."""
    m, n = cost.shape
    model = gp.Model("gap")
    model.Params.OutputFlag = 0
    x = model.addVars(m, n, vtype=GRB.BINARY, name="x")
    model.setObjective(
        gp.quicksum(cost[i, j] * x[i, j] for i in range(m) for j in range(n)),
        GRB.MINIMIZE,
    )
    add_job_assignment_constraints(model, x, m, n)
    add_capacity_constraints(model, x, a, b)
    return model, x

def solve_gap_mip(cost: np.ndarray, a: np.ndarray, b: np.ndarray,
                  time_limit: float = 60.0) -> tuple[np.ndarray, float, float]:
    """Solve GAP to optimality or time limit.

    Returns (assign, objective, mip_gap) where assign[j] is job j's agent.
    """
    m, n = cost.shape
    model, x = build_gap_model(cost, a, b)
    model.Params.TimeLimit = time_limit
    model.Params.MIPGap = 1e-6
    model.optimize()
    solved = model.Status == GRB.OPTIMAL or (
        model.Status == GRB.TIME_LIMIT and model.SolCount > 0
    )
    if not solved:
        raise RuntimeError(f"no solution: status {model.Status}")
    assign = np.array(
        [max(range(m), key=lambda i: x[i, j].X) for j in range(n)], dtype=int
    )
    return assign, float(model.ObjVal), float(model.MIPGap)

cost = np.array([[8.0, 6.0, 5.0, 7.0],
                 [6.0, 7.0, 8.0, 5.0]])
a = np.array([[3, 2, 4, 3],
              [3, 3, 3, 2]])
b = np.array([7, 6])
assign, obj, gap = solve_gap_mip(cost, a, b)
print(assign, obj, gap)
# Expected: assign [1 0 0 1], objective 22.0, gap 0.0
# (every job at its cheapest agent happens to fit: loads 6  tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Generate a min-cost GAP instance (cost, a, b), Chu-Beasley classes.

    'c': independent uniform cost and resource.
    'd': cost = 111 - a + noise (inverse correlation, hard).
    'e': resource log-distributed, cost ~ 1000/a (hardest).
    Capacities are 80% of average load; feasibility is NOT guaranteed by
    construction (GAP feasibility itself is NP-complete) - let the solver
    or validator detect infeasible draws.
    """
    rng = np.random.default_rng(seed)
    if gap_class == "c":
        a = rng.integers(5, 26, size=(n_agents, n_jobs))
        cost = rng.integers(10, 51, size=(n_agents, n_jobs))
    elif gap_class == "d":
        a = rng.integers(1, 101, size=(n_agents, n_jobs))
        cost = 111 - a + rng.integers(-10, 11, size=(n_agents, n_jobs))
    elif gap_class == "e":
        u = rng.random((n_agents, n_jobs))
        a = np.maximum(1, (1.0 - 10.0 * np.log(u)).astype(int))
        cost = np.maximum(1, (1000.0 / a - 10.0 * rng.random((n_agents, n_jobs))).astype(int))
    else:
        raise ValueError(f"unknown gap_class {gap_class!r}")
    b = np.maximum(
        np.ceil(0.8 * a.sum(axis=1) / n_agents).astype(int),
        a.max(axis=1),                       # each job must fit somewhere alone
    )
    return cost.astype(float), a.astype(int), b.astype(int)

cost, a, b = generate_gap_instance(3, 12, "c", seed=42)
print(cost.shape, a.shape, b)
# Expected: (3, 12) (3, 12) and a 3-entry capacity vector near
# 0.8 * (row resource sum) / 3, e.g. around 45-55 for these parameters

Independent feasibility and objective validator

Never report a solution straight from a solver or heuristic. Recompute feasibility and objective from the raw data with code that shares nothing with the model. This catches indexing bugs, sense errors, and stale data.

import numpy as np

def validate_lap_solution(col_of_row: np.ndarray, cost: np.ndarray) -> float:
    """Check a square-LAP solution is a permutation; return recomputed cost."""
    n = cost.shape[0]
    if sorted(col_of_row.tolist()) != list(range(n)):
        raise ValueError("not a permutation: a column is reused or missing")
    return float(cost[np.arange(n), col_of_row].sum())

def validate_gap_solution(assign: np.ndarray, cost: np.ndarray, a: np.ndarray,
                          b: np.ndarray) -> tuple[bool, float, list[str]]:
    """Independent feasibility + objective check for a GAP solution.

    assign[j] is the agent of job j. Recomputes loads and cost from raw
    data only. Returns (feasible, objective, issue messages).
    """
    m, n = cost.shape
    issues: list[str] = []
    if assign.shape != (n,):
        return False, float("nan"), [f"assign shape {assign.shape} != ({n},)"]
    if ((assign = m)).any():
        return False, float("nan"), ["agent index out of range"]
    jobs = np.arange(n)
    load = np.bincount(assign, weights=a[assign, jobs], minlength=m)
    for i in np.flatnonzero(load > b + 1e-9):
        issues.append(f"agent {i}: load {load[i]:.0f} > capacity {b[i]:.0f}")
    objective = float(cost[assign, jobs].sum())
    return not issues, objective, issues

cost = np.array([[8.0, 6.0, 5.0, 7.0],
                 [6.0, 7.0, 8.0, 5.0]])
a = np.array([[3, 2, 4, 3], [3, 3, 3, 2]])
b = np.array([7, 6])
print(validate_gap_solution(np.array([1, 0, 0, 1]), cost, a, b))
# Expected: (True, 22.0, [])
print(validate_gap_solution(np.array([0, 0, 0, 0]), cost, a, b))
# Expected: (False, 26.0, ['agent 0: load 12 > capacity 7'])

Generalized Assignment: Lagrangian Bound and Local Search

Lagrangian relaxation with knapsack subproblems

Dualiz

…

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: hajibabaie
• Source: hajibabaie/combinatorial-optimization-skills
• License: MIT

Install and usage instructions live in the source repository linked above.