Set Covering Packing Partitioning

Set Covering, Packing, and Partitioning

You are an expert in integer programming and combinatorial optimization with deep experience in covering-type models. This skill covers the set covering problem (SCP), set packing, and set partitioning: exact MIP formulations, the greedy ln(n)-approximation, LP rounding, Lagrangian-based heuristics in the Caprara-Fischetti-Toth tradition, GRASP, and the crew-scheduling and column-generation context where these models dominate. Use the framework below to classify the variant first, then pick a solution path that matches instance size and optimality needs.

Initial Assessment

Establish these facts before writing any model or code:

• Row semantics. Must each row be hit at least once (covering), exactly once (partitioning), or at most once (packing)? This single choice changes feasibility behavior, LP strength, and the right algorithm more than anything else.
• Direction of the objective. Covering and partitioning minimize cost; packing maximizes profit. Confirm the user is not mixing the two (e.g., "maximize coverage" is a different model: maximal covering location, a budgeted variant).
• Instance size and density. How many rows m, columns n, and what fraction of entries are nonzero? OR-Library SCP instances run from 200 x 1000 at 2% density to unicost instances at 5,000+ columns. Density drives both memory layout and LP difficulty.
• Are the columns explicit or implicit? If columns are routes, crew pairings, or patterns that must be generated, this is a column-generation master problem, not a static MIP. Route the conversation accordingly.
• Cost structure. Unicost (all c_j = 1) instances behave very differently from general-cost ones: LP bounds are weaker and greedy is less informative. Ask.
• Side constraints. Base capacities in crew scheduling, cardinality limits, mutual-exclusion pairs. These break pure SCP structure and may push toward a plain MIP or branch-and-price.
• Over-covering tolerance. In practice, can a row be covered twice (a deadheading crew, a doubly visited customer)? If yes, replace fragile partitioning constraints with covering plus a penalty.
• Solver availability and time budget. Gurobi licensed? Seconds or hours? Proven optimal or "good cover fast"? Gurobi solves most general-cost SCPs of moderate size outright; large unicost and partitioning instances are where heuristics earn their keep.
• Validation requirements. Agree up front that every reported solution passes an independent feasibility and objective check, separate from the model code.

Problem Definitions and Formulations

Ground set $I = {1, \dots, m}$ of rows (elements, tasks, flight legs). A family of $n$ columns (subsets) $S_j \subseteq I$, each with cost or profit $c_j > 0$. The incidence matrix $A \in {0,1}^{m \times n}$ has $a_{ij} = 1$ iff $i \in S_j$. Binary variable $x_j = 1$ iff column $j$ is selected.

Set covering (SCP) — every row hit at least once, minimize cost:

$$ \min ; \sum_{j=1}^{n} c_j x_j \quad \text{s.t.} \quad \sum_{j:, i \in S_j} x_j \ge 1 ;; \forall i \in I, \qquad x \in {0,1}^n $$

Set packing — rows hit at most once, maximize profit:

$$ \max ; \sum_{j=1}^{n} c_j x_j \quad \text{s.t.} \quad \sum_{j:, i \in S_j} x_j \le 1 ;; \forall i \in I, \qquad x \in {0,1}^n $$

Set partitioning — every row hit exactly once, minimize cost:

$$ \min ; \sum_{j=1}^{n} c_j x_j \quad \text{s.t.} \quad \sum_{j:, i \in S_j} x_j = 1 ;; \forall i \in I, \qquad x \in {0,1}^n $$

All three are NP-hard; SCP and exact cover (partitioning feasibility) are among Karp's original 21 problems (Karp 1972, "Reducibility Among Combinatorial Problems").

Variant selection

Relationships worth exploiting:

• A partition is exactly a solution that is simultaneously a cover and a packing.
• If over-covering is physically harmless (free disposal), solve the covering relaxation of a partitioning model: every cover can be interpreted as a partition plus tolerated duplicates, and the covering LP is far better behaved (less degenerate, always feasible when every row has a column).
• Set packing with rows = edges and columns = vertices is maximum-weight independent set, so packing inherits its inapproximability: no $n^{1-\varepsilon}$ approximation unless P = NP (Hastad 1999, "Clique is hard to approximate"; Zuckerman 2007 derandomization).
• Partitioning feasibility alone is NP-complete (exact cover), so heuristics need an explicit feasibility strategy, not just an objective-improvement loop.

Approximability and bound quality

Decision guidance:

• Use the MIP directly when n is up to a few hundred thousand explicit columns and costs are general. Modern solvers exploit the structure well.
• Use greedy / GRASP / Lagrangian heuristics for unicost instances, very large n, tight time budgets, or as upper bounds and MIP starts.
• Use column generation / branch-and-price when columns are implicit (pairings, routes, patterns); see the Advanced Techniques section.
• Reformulate partitioning as penalized covering whenever the application tolerates over-coverage; this is standard practice in crew scheduling.

Exact Models in Gurobi

One builder function per constraint family keeps the model auditable: each family can be unit-tested on a tiny instance, constraint names make IIS output readable when a partitioning model turns out infeasible, and switching variant means swapping one builder. For partitioning instances where finding any feasible point is hard, set MIPFocus = 1 and consider seeding a cover-based start (see Practical Challenges).

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


def add_cover_constraints(model: gp.Model, x: gp.tupledict, A: np.ndarray) -> None:
    """Covering rows: each row is hit at least once."""
    for i in range(A.shape[0]):
        cols = np.flatnonzero(A[i])
        model.addConstr(gp.quicksum(x[j] for j in cols) >= 1, name=f"cover[{i}]")


def add_packing_constraints(model: gp.Model, x: gp.tupledict, A: np.ndarray) -> None:
    """Packing rows: each row is hit at most once."""
    for i in range(A.shape[0]):
        cols = np.flatnonzero(A[i])
        model.addConstr(gp.quicksum(x[j] for j in cols) <= 1, name=f"pack[{i}]")


def add_partition_constraints(model: gp.Model, x: gp.tupledict, A: np.ndarray) -> None:
    """Partitioning rows: each row is hit exactly once."""
    for i in range(A.shape[0]):
        cols = np.flatnonzero(A[i])
        model.addConstr(gp.quicksum(x[j] for j in cols) == 1, name=f"part[{i}]")


def build_set_system_model(
    A: np.ndarray, c: np.ndarray, mode: str
) -> tuple[gp.Model, gp.tupledict]:
    """Build the SCP / set packing / set partitioning MIP over incidence matrix A."""
    n = A.shape[1]
    model = gp.Model(f"set-{mode}")
    model.Params.OutputFlag = 0
    x = model.addVars(n, vtype=GRB.BINARY, name="x")
    sense = GRB.MAXIMIZE if mode == "packing" else GRB.MINIMIZE
    model.setObjective(gp.quicksum(float(c[j]) * x[j] for j in range(n)), sense)
    builders = {
        "cover": add_cover_constraints,
        "packing": add_packing_constraints,
        "partition": add_partition_constraints,
    }
    builders[mode](model, x, A)
    return model, x


def solve_set_system(
    A: np.ndarray, c: np.ndarray, mode: str, time_limit: float = 60.0
) -> dict[str, object]:
    """Solve the chosen variant and extract the incumbent; raise if none exists."""
    model, x = build_set_system_model(A, c, mode)
    model.Params.TimeLimit = time_limit
    model.Params.MIPGap = 1e-6
    model.optimize()
    usable = model.Status == GRB.OPTIMAL or (
        model.Status == GRB.TIME_LIMIT and model.SolCount > 0
    )
    if not usable:
        raise RuntimeError(f"no solution available, Gurobi status {model.Status}")
    sol = np.array([x[j].X > 0.5 for j in range(A.shape[1])])
    return {"x": sol, "objective": model.ObjVal, "bound": model.ObjBound,
            "gap": model.MIPGap, "status": model.Status}


# Tiny instance: 4 rows, 5 columns S1..S5 with S3 = {1,3}, S4 = {2,4}.
A = np.array([[1, 0, 1, 0, 1], [1, 0, 0, 1, 1],
              [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]], dtype=bool)
c = np.array([3.0, 3.0, 2.0, 2.0, 5.0])
for mode in ("cover", "partition", "packing"):
    res = solve_set_system(A, c, mode)
    print(mode, res["objective"], sorted(np.flatnonzero(res["x"])))
# Expected: cover 4.0 [2, 3] / partition 4.0 [2, 3] / packing 6.0 [0, 1]

The same data yields three different optima: columns S3 and S4 form both the cheapest cover and the cheapest partition (cost 4), while the most profitable packing takes S1 and S2 (profit 6). Keeping the three builders separate makes such comparisons a two-line experiment.

Instance Generation and Independent Validation

Benchmark conventions follow Beasley (1987, 1990; OR-Library): classes scp4x-scp6x and scpa-scpe have densities of 2-20% and costs uniform in [1, 100]; the scpnr and unicost classes (including Steiner triple systems) are markedly harder for both LP bounds and heuristics. The generator below reproduces that style with two guarantees: every row is coverable by at least two columns (so covering is feasible and no column is trivially forced), and no column is empty.

import numpy as np


def generate_scp_instance(
    m: int, n: int, density: float = 0.05, cost_mode: str = "random", seed: int = 0
) -> tuple[np.ndarray, np.ndarray]:
    """Random set-system instance in OR-Library style (Beasley 1987).

    cost_mode: 'unicost' (all 1), 'random' (uniform integers in [1, 100]),
    or 'size_correlated' (cost grows with column size, which weakens greedy).
    """
    rng = np.random.default_rng(seed)
    A = rng.random((m, n)) < density
    # Repair empty columns: give each one a random row.
    empty_cols = np.flatnonzero(~A.any(axis=0))
    A[rng.integers(0, m, size=empty_cols.size), empty_cols] = True
    # Repair rows covered fewer than twice.
    for i in np.flatnonzero(A.sum(axis=1) < 2):
        need = 2 - int(A[i].sum())
        candidates = np.flatnonzero(~A[i])
        A[i, rng.choice(candidates, size=need, replace=False)] = True
    if cost_mode == "unicost":
        c = np.ones(n)
    elif cost_mode == "random":
        c = rng.integers(1, 101, size=n).astype(float)
    else:  # size_correlated
        c = np.ceil(A.sum(axis=0) * rng.uniform(0.8, 1.2, size=n) * 10.0)
    return A, c


A, c = generate_scp_instance(m=50, n=200, density=0.05, cost_mode="random", seed=42)
print(A.shape, int(A.sum(axis=1).min()), float(c.min()), float(c.max()))
# Expected: (50, 200) with minimum row coverage >= 2 and costs inside [1, 100].

The validator recomputes feasibility and objective from the raw incidence matrix only. It shares no code with the model builders, so a bug in constraint construction cannot hide itself. Run it on every solution you report, from any method.

import numpy as np


def validate_solution(
    A: np.ndarray, c: np.ndarray, x: np.ndarray, mode: str
) -> dict[str, object]:
    """Independent feasibility + objective check from raw data; trusts no solver."""
    x = np.asarray(x, dtype=bool)
    coverage = A[:, x].sum(axis=1)
    if mode == "cover":
        violated = np.flatnonzero(coverage < 1)
    elif mode == "partition":
        violated = np.flatnonzero(coverage != 1)
    elif mode == "packing":
        violated = np.flatnonzero(coverage > 1)
    else:
        raise ValueError(f"unknown mode {mode!r}")
    return {
        "feasible": violated.size == 0,
        "objective": float(c[x].sum()),
        "violated_rows": violated.tolist(),
        "n_selected": int(x.sum()),
        "max_coverage": int(coverage.max()) if coverage.size else 0,
    }


A = np.array([[1, 0, 1, 0, 1], [1, 0, 0, 1, 1],
              [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]], dtype=bool)
c = np.array([3.0, 3.0, 2.0, 2.0, 5.0])
report = validate_solution(A, c, np.array([0, 0, 1, 1, 0]), "partition")
print(report["feasible"], report["objective"], report["violated_rows"])
# Expected: True 4.0 [] -- columns S3 and S4 partition all four rows at cost 4.

The max_coverage field is a cheap diagnostic: a covering solution with max_coverage well above 1 usually still contains redundant columns and should go through reverse delete before being reported.

Approximation Algorithms: Greedy and LP Rounding

The greedy heuristic (Chvatal 1979) repeatedly selects the column minimizing $c_j / |\text{newly covered rows}|$. Its cost is at most $H(d) \cdot \mathrm{OPT}$ with $d = \max_j |S_j|$ and $H(d) \le \ln d + 1$, and this is essentially the best possible for polynomial algorithms: Feige (1998) showed no $(1-o(1))\ln n$ approximation exists unless P = NP. Greedy output often contains redundant columns picked early, so always follow it with reverse delete (drop expensive columns whose rows remain covered).

import numpy as np


def greedy_set_cover(A: np.ndarray, c: np.ndarray) -> np.ndarray:
    """Chvatal (1979) greedy: best cost per newly covered row, until all covered.

    Guarantee: cost <= H(d) * OPT with d = largest column size. Assumes every
    row is covered by at least one column.
    """
    m, n = A.shape
    uncovered = np.ones(m, dtype=bool)
    selected = np.zeros(n, dtype=bool)
    while uncovered.any():
        gain = A[uncovered].sum(axis=0).astype(float)
        gain[selected] = 0.0
        ratio = np.full(n, np.inf)
        pos = gain > 0
        ratio[pos] = c[pos] / gain[pos]
        j = int(np.argmin(ratio))
        selected[j] = True
        uncovered &= ~A[:, j]
    return selected


def remove_redundant_columns(
    A: np.ndarray, c: np.ndarray, selected: np.ndarray
) -> np.ndarray:
    """Reverse delete: drop the most expensive redundant columns first."""
    selected = selected.copy()
    cover = A[:, selected].sum(axis=1)
    for j in sorted(np.flatnonzero(selected), key=lambda k: -c[k]):
        rows = A[:, j]
        if np.all(cover[rows] >= 2):
            selected[j] = False
            cover[rows] -= 1
    return selected


A = np.array([[1, 0, 1, 0, 1], [1, 0, 0, 1, 1],
              [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]], dtype=bool)
c = np.array([3.0, 3.0, 2.0, 2.0, 5.0])
x = remove_redundant_columns(A, c, greedy_set_cover(A, c))
print(sorted(np.flatnonzero(x)), float(c[x].sum()))
# Expected: [2, 3] 4.0 -- the MIP optimum; the a-priori guarantee was H(4) ~ 2.08 * OPT.

LP rounding (Hochbaum 1982) gives a complementary guarantee tied to row frequency rather than column size. Solve the LP relaxation, then select every column with $x_j^{LP} \ge 1/f$, where $f$ is the maximum number of columns through any row. Each covering row has at most $f$ fractional variables summing to at least 1, so at least one of them clears the threshold: the rounded solution is always feasible, and its cost is at most $f \cdot \mathrm{OPT}$. With $f = 2$ this is exactly the classic vertex-cover 2-approximation. The LP value is also the standard lower bound to report alongside any heuristic. Randomized rounding (select column $j$ with probability $\min(1, \alpha x_j^{LP})$, $\alpha = O(\ln m)$, repeat until feasible) usually does better in cost but needs the repetition loop.

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


def lp_rounding_cover(A: np.ndarray, c: np.ndarray) -> tuple[np.ndarray, float]:
    """Hochbaum (1982) deterministic LP rounding: an f-approximation for SCP.

    f = max number of columns covering any single row. Returns (solution, LP bound).
    """
    m, n = A.shape
    model = gp.Model("scp-lp")
    model.Params.OutputFlag = 0
    x = model.addVars(n, lb=0.0, ub=1.0, name="x")
    for i in range(m):
        cols = np.flatnonzero(A[i])
        model.addConstr(gp.quicksum(x[j] for j in cols) >= 1, name=f"cover[{i}]")
    model.setObjective(gp.quicksum(float(c[j]) * x[j] for j in range(n)), GRB.MINIMIZE)
    model.optimize()
    if model.Status != GRB.OPTIMAL:
        raise RuntimeError(f"LP not solved to optimality, status {model.Status}")
    x_lp = np.array([x[j].X for j in range(n)])
    f = int(A.sum(axis=1).max())
    rounded = x_lp >= 1.0 / f - 1e-9
    return rounded, float(model.ObjVal)


A = np.array([[1, 0, 1, 0, 1], [1, 0, 0, 1, 1],
              [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]], dtype=bool)
c = np.array([3.0, 3.0, 2.0, 2.0, 5.0])
x_round, lp_bound = lp_rounding_cover(A, c)
print(sorted(np.flatnonzero(x_round)), float(c[x_round].sum()), lp_bound)
# Expected: [2, 3] 4.0 4.0 -- the LP is integral here, so rounding returns the optimum.

Practical reading of the two guarantees: greedy wins when columns are large and rows are in many columns (big $f$, moderate $d$); LP rounding wins when each row appears in few columns (small $f$), as in vertex-cover-like structures. On general-cost OR-Library instances, both typically land within a few percent of the LP bound, and the LP bound itself is within 1-2% of the integer optimum.

GRASP and the Caprara-Fischetti-Toth Lagrangian Heuristic

GRASP is a natural fit for SCP — it was introduced on this very problem (Feo & Resende 1989, "A probabilistic heuristic for a computationally difficult set covering problem"). Each iteration builds a cover with a restricted candidate list (RCL) over the greedy cost-effectiveness ratio, then cleans it with reverse delete. The implementation below stays under 80 lines; the grasp skill covers reactive alpha tuning, RCL design choices, and path relinking between elite covers, all of which drop into this loop unchanged.

import numpy as np


def _rcl_construction(
    A: np.ndarray, c: np.ndarray, alpha: float, rng: np.random.Generator
) -> np.ndarray:
    """Greedy randomized construction: sample uniformly from the RCL of columns
    whose cost-effectiveness is within alpha of the best."""
    m, n = A.shape
    uncovered = np.ones(m, dtype=bool)
    selected = np.zeros(n, dtype=bool)
    while uncovered.any():
        gain = A[uncovered].sum(axis=0).astype(float)
        gain[selected] = 0.0
        ratio = np.full(n, np.inf)
        pos = gain > 0
        ratio[pos] = c[pos] / gain[pos]
        r_min, r_max = ratio[pos].min(), ratio[pos].max()
        threshold = r_min + alpha * (r_max - r_min)
        rcl = np.flatnonzero(ratio <= threshold + 1e-12)
        j = int(rng.choice(rcl))
        selected[j] = True
        uncovered &= ~A[:, j]
    return selected


def _drop_redundant(A: np.ndarray, c: np.ndarray, selected: np.ndarray) -> np.ndarray:
    """Reverse delete: drop selected columns (most expensive first) not needed."""
    selected = selected.copy()
    cover = A[:, selected].sum(axis=1)
    for j in sorted(np.flatnonzero(selected), key=lambda k: -c[k]):
        rows = A[:, j]
        if np.all(cover[rows] >= 2):
            selected[j] = False
            cover[rows] -= 1
    return selected


def grasp_set_cover(
    A: np.ndarray, c: np.ndarray, n_iters: int = 200, alpha: float = 0.3, seed: int = 0
) -> tuple[np.ndarray, float]:
    """GRASP for SCP (Feo & Resende 1989): RCL construction + redundancy removal.

    alpha = 0 is pure greedy; alpha = 1 is uniformly random construction.
    See the grasp skill for reactive alpha and path relinking extensions.
    """
    rng = np.random.default_rng(seed)
    best_x = np.zeros(A.shape[1], dtype=bool)
    best_cost = np.inf
    for _ in range(n_iters):
        x = _drop_redundant(A, c, _rcl_construction(A, c, alpha, rng))
        cost = float(c[x].sum())
        if cost < best_cost:
            best_x, best_cost = x.copy(), cost
    return best_x, best_cost


A = np.array([[1, 0, 1, 0, 1], [1, 0, 0, 1, 1],
              [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]], dtype=bool)
c = np.array([3.0, 3.0, 2.0, 2.0, 5.0])
x_best, cost_best = grasp_set_cover(A, c, n_iters=50, alpha=0.3, seed=7)
print(sorted(np.flatnonzero(x_best)), cost_best)
# Expected: [2, 3] 4.0 -- GRASP finds the optimal cover {S3, S4} on this instance.

The strongest classical SCP heuristic is Lagrangian-based: Caprara, Fischetti & Toth (1999), "A heuristic method for the set covering problem", won the FASTER competition on instances with millions of columns. Dualize the covering constraints with multipliers $u \ge 0$:

$$ L(u) = \sum_{i} u_i + \sum_{j} \min\bigl(0,; c_j - \textstyle\sum_{i \in S_j} u_i\bigr), \qquad x_j(u) = 1 \iff c_j - \textstyle\sum_{i \in S_j} u_i < 0 . $$

Maximize $L(u)$ by subgradient ascent ($g_i = 1 - \sum_j a_{ij} x_j(u)$, step $\lambda (UB - L(u)) / \lVert g \rVert^2$), and periodically recover primal covers by running greedy on the Lagrangian reduced costs — near-optimal multipliers make reduced costs an excellent column-quality signal. Because the Lagrangian subproblem has the integrality property, the best achievable bound equals the LP bound, but it arrives orders of magnitude faster than solving the LP on huge instances. The lagrangian-relaxation skill covers step-size rules and convergence in depth.

import numpy as np


def _greedy_with_costs(A: np.ndarray, score: np.ndarray) -> np.ndarray:
    """Greedy cover driven by an arbitrary nonnegative score vector as costs."""
    m, n = A.shape
    uncovered = np.ones(m, dtype=bool)
    selected = np.zeros(n, dtype=bool)
    while uncovered.any():
        gain = A[uncovered].sum(axis=0).astype(float)
        gain[selected] = 0.0
        ratio = np.where(gain > 0, score / np.maximum(gain, 1e-12), np.inf)
        j = int(np.argmin(ratio))
        selected[j] = True
        uncovered &= ~A[:, j]
    return selected


def lagrangian_scp(
    A: np.ndarray, c: np.ndarray, n_iters: int = 300, seed: int = 0
) -> dict[str, float]:
    """CFT-style subgradient + Lagrangian greedy (Caprara, Fischetti & Toth 1999).

    Returns the best lower bound, best validated upper bound, and gap percent.
    """
    m, n = A.shape
    Af = A.astype(float)
    cost_per_row = c / A.sum(axis=0)          # cost per covered row, per column
    u = np.array([cost_per_row[A[i]].min() for i in range(m)])
    x0 = _greedy_with_costs(A, c.astype(float))
    best_x, best_ub = x0, float(c[x0].sum())
    best_lb, lam, stall = 0.0, 2.0, 0
    for it in range(n_iters):
        reduced = c - Af.T @ u
        x_u = reduced < 0.0
        lb = float(u.sum() + reduced[x_u].sum())
        if lb > best_lb + 1e-9:
            best_lb, stall = lb, 0
        else:
            stall += 1
        if stall == 20:
            lam, stall = lam / 2.0, 0
        g = 1.0 - Af @ x_u
        sq = float(g @ g)
        if sq < 1e-12:
            break                              # x(u) feasible: multipliers optimal
        if it % 10 == 0:                       # primal recovery on reduced costs
            x_h = _greedy_with_costs(A, np.maximum(reduced, 0.0) + 1e-6)
            ub = float(c[x_h].sum())
            if ub < best_ub:
                best_ub, best_x = ub, x_h
        u = np.maximum(0.0, u + lam * (best_ub - lb) / sq * g)
    gap = 100.0 * (best_ub - best_lb) / best_ub
    return {"best_lb": best_lb, "best_ub": best_ub,
            "gap_pct": gap, "n_selected": float(best_x.sum())}


A = np.array([[1, 0, 1, 0, 1], [1, 0, 0, 1, 1],
              [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]], dtype=bool)
c = np.array([3.0, 3.0, 2.0, 2.0, 5.0])
res = lagrangian_scp(A, c, n_iters=100, seed=1)
print(round(res["best_lb"], 2), res["best_ub"], round(res["gap_pct"], 2))
# Expected: 4.0 4.0 0.0 -- the Lagrangian dual matches the LP bound and closes the gap.

Advanced Techniques

Lagrangian cores and column fixing

The full CFT method wraps the subgradient loop above in two scaling devices. First, a core problem: keep only the 5 or so columns with the smallest Lagrangian reduced cost per row (plus columns in the incumbent), run everything on this core, and periodically re-price the full column set against the current multipliers to refresh it. This is what makes million-column instances tractable in memory and time. Second, reduced-cost fixing: with lower bound $LB$ and incumbent $UB$, any column with $c_j - \sum_{i \in S_j} u_i > UB - LB$ cannot appear in an improving solution and is fixed to zero; conversely, CFT fixes to one the columns that persist across several near-best Lagrangian greedy solutions, then re-optimizes the residual problem. Both devices transfer directly to any Lagrangian heuristic you build on covering-type models.

Crew scheduling as set partitioning

Airline and rail crew pairing is the canonical industrial set partitioning problem (Garfinkel & Nemhauser 1969; Barnhart et al. 1998, "Branch-and-Price: Column Generation for Solving Huge Integer Programs"). Rows are flight legs or trips in a planning horizon; a column is a legal pairing — a multi-day sequence of legs respecting duty-time, rest, and connection rules — with cost equal to guaranteed pay minus flying credit. The model is $\min c^\top x$, $Ax = \mathbf{1}$, plus side constraints such as crew-base capacities $\sum_{j \in B_b} x_j \le \kappa_b$. Practitioners almost always solve the covering relaxation $Ax \ge \mathbf{1}$ with an over-cover penalty: a doubly covered leg means one crew deadheads on it, which is operationally legal and priced at a real cost. This single modeling choice removes most partitioning infeasibility and LP degeneracy pain. Calibrate the penalty to the true deadhead cost, not to an arbitrary big-M, or the optimizer will trade real flying for fake savings.

The column generation connection

In crew pairing, vehicle routing, and cutting stock, the column set is implicit and astronomically large, so the partitioning/covering model is solved as a restricted master problem: keep a manageable subset of columns, get duals $u_i$ from the LP, and price new columns by solving $\min_j c_j - \sum_{i \in S_j} u_i$ over the implicit set — a resource-constrained shortest path for pairings and routes, a knapsack for cutting patterns. Negative reduced cost columns enter; repeat to LP optimality, then branch with rules that preserve pricing structure (Ryan & Foster 1981: branch on a row pair (s, t) being covered by the same or different columns). The Lagrangian multipliers from the CFT loop are interchangeable with LP duals here, which is why the two literatures cross-reference constantly. The column-generation skill covers the master-pricing loop, stabilization, and branch-and-price mechanics; the vehicle-routing-problem skill shows the route-based master in full.

Preprocessing and problem reductions

Apply these dominance rules before any algorithm; on OR-Library instances they often remove 10-50% of the matrix (Beasley 1987):

• Row domination. If the column set covering row $k$ is a subset of the columns covering row $i$, any solution covering $k$ also covers $i$: delete row $i$.
• Column domination (covering only). If $S_{j'} \supseteq S_j$ and $c_{j'} \le c_j$, delete column $j$. Never apply this to partitioning — the smaller column may be needed for exactness.
• Forced columns. A row covered by a single column fixes that column to 1; remove all rows it covers (and, in partitioning, delete every column that intersects them — this cascade is powerful).
• Duplicate columns. Keep the cheapest of each identical column; in packing, keep the most profitable.

Iterate the rules to a fixed point: each fixing can enable new dominations.

Practical Challenges

The set partitioning model is infeasible, or the solver finds no incumbent for hours. Partitioning feasibility is itself NP-complete, so do not expect the solver to stumble into it. Add one artificial slack column per row with a cost above any plausible over-cover penalty, or switch to the covering relaxation with calibrated penalties. Seed a MIP start from a greedy cover with duplicates removed. Set MIPFocus = 1 to bias the search toward feasibility.

The partitioning LP is massively degenerate and dual values swing wildly between iterations. Equality rows with many identical-cost columns create degenerate vertices. Use barrier without crossover (Method = 2, Crossover = 0) when you only need dual prices for pricing or fixing, perturb costs by tiny amounts, or stabilize duals around a Lagrangian center as in CFT.

Unicost instances give weak LP bounds and indifferent greedy choices. With all $c_j = 1$ the cost-effectiveness ratio collapses to coverage counts and ties are everywhere. Break ties randomly across restarts (GRASP does this for free), use Lagrangian multipliers to differentiate columns, and report gaps against the Lagrangian bound rather than the raw LP when the LP is too slow.

The incidence matrix does not fit in dense memory. A 10,000 x 1,000,000 instance is 10 GB dense but a few hundred MB sparse. Store columns as index lists or scipy.sparse CSC; build Gurobi constraints from the sparse structure directly. The dense-numpy code in this skill is for clarity at moderate sizes — swap the matrix products for sparse ones without changing the algorithms.

Branching on fractional x_j stalls branch-and-bound on partitioning models. Fixing a single column variable to 0 barely changes the LP, producing huge, unbalanced trees. Use Ryan-Foster row-pair branching (same/different), which splits the solution space evenly and remains compatible with column generation pricing.

Heuristic covers carry hidden redundancy. Greedy and GRASP constructions routinely include early columns that later picks make unnecessary; reporting them overstates cost by several percent. Always run reverse delete, and check the validator's max_coverage output — values above 2 on most rows signal a sloppy construction loop.

Over-cover penalties distort the objective. If the penalty for double-covering a row is far above its real operational cost, the model buys expensive exact partitions it does not need; far below, it spams duplicates. Price the penalty at the measured cost of the real-world event (deadhead, double visit) and run a sensitivity sweep over a 2x range to confirm solution stability.

Identical and dominated columns inflate the model and create symmetry. Column generators (route enumerators, pairing builders) frequently emit duplicates and dominated variants, multiplying symmetric optima and slowing the MIP. Run the preprocessing reductions above after every generation batch, not just once.

Tools & Libraries

Output Format

A complete answer to a covering/packing/partitioning task contains:

1. Variant statement. One sentence fixing the constraint sense (>= / = / <=), the objective direction, and any side constraints, plus the m x n size and density of the instance.
2. Method and bound. What was run (MIP, greedy, LP rounding, GRASP, Lagrangian) with parameters and seed, the best solution cost, and the best available lower (upper, for packing) bound with the resulting gap.
3. Validator confirmation. The output of validate_solution on the reported solution — never report a solution that has not passed it.
4. Solution artifact. The selected column indices (and, where meaningful, the row-to-column assignment), written to a file for nontrivial sizes.

A reporting template:

Instance      : scp45 (m=200, n=1000, density=2.0%, costs in [1,100])
Variant       : set covering (Ax >= 1, minimize)
Method        : CFT Lagrangian heuristic, 300 subgradient iters, seed 42
Best cost     : 512.0       (validator: feasible=True, objective=512.0)
Lower bound   : 508.3       (Lagrangian dual, = LP bound up to convergence)
Gap           : 0.73 %
Time          : 1.8 s
Artifact      : results/scp45_seed42_solution.json (87 columns selected)

For method comparisons, add a table with one row per (method, seed) and columns {cost, bound, gap %, time, n_selected}, aggregated over seeds with mean and best. State explicitly which preprocessing reductions were applied, since they change the effective instance.

Questions to Ask

• Must every row be hit at least once, exactly once, or at most once — and what does an over-covered row mean physically in your application?
• Are the columns given explicitly as data, or generated from rules (routes, pairings, patterns)? Roughly how many are there?
• What are m, n, and the matrix density? Unicost or general costs?
• Are there side constraints beyond the row constraints (capacities, cardinality, incompatible column pairs)?
• Do you need a proven optimum, a bounded gap, or just a good solution within a time budget — and what is that budget?
• Is a Gurobi license available, or should the solution path use open-source tools?
• Do you have benchmark or historical instances to calibrate against, or should we generate synthetic ones with controlled density and cost structure?
• How will solutions be consumed downstream — do you need the row-to-column assignment, or just the selected columns and cost?

Related Skills

• grasp — when the greedy-randomized construction here needs reactive alpha tuning, RCL design depth, or path relinking between elite covers.
• column-generation — when columns are implicit and the covering/partitioning model is the restricted master of a pricing loop or branch-and-price.
• milp-modeling-gurobi — when the exact models need richer gurobipy patterns: parameters, status handling, warm starts, and constraint-builder organization.
• lagrangian-relaxation — when the CFT-style subgradient loop needs step-size rules, convergence diagnostics, or different dualization choices.
• vehicle-routing-problem — when the set-partitioning master's columns are vehicle routes and the application is routing rather than crews or patterns.