Big-M constants come from declared bounds#

declared_bounds_bigm.ipynb Open In Colab Open In Deepnote Open In Kaggle Open In Gradient Open In SageMaker Studio Lab Powered by AMPL

Description: How MP sizes the big-M constants of logical constraints from declared variable bounds, why naming an expression and declaring its bound tightens the LP relaxation, and how to check this with the root bound and the exported solver model

Tags: modeling-tips, mp, big-m, reinforced-formulations, formulation-comparison, facility-location, milp, amplpy, highs

Notebook author: Gleb Belov <gleb@ampl.com>

References:

  1. MP modeling guide, Efficiency considerations and Numerical accuracy: https://mp.ampl.com/modeling-efficiency.html, https://mp.ampl.com/modeling-numerics.html

  2. Vielma, J.P. (2015). Mixed integer linear programming formulation techniques. SIAM Review 57(1), 3-57.

  3. Bonami, P., Lodi, A., Tramontani, A., Wiese, S. (2015). On mathematical programming with indicator constraints. Mathematical Programming 151(1), 191-223.

# Install dependencies
%pip install -q amplpy numpy pandas

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

SOLVER_MILP = (
    "highs"  # any MP-based MILP solver: highs, gurobi, xpress, cplex, copt, scip, ...
)
# HiGHS has no native indicator constraints, so MP linearizes every implication
# with a big-M constant -- which is what this notebook inspects. A solver with
# native indicators may keep some implications as indicators instead.

The point#

When a solver has no native support for a logical constraint such as

\[\text{Open}_i = 0 \;\Longrightarrow\; \sum_j \text{Ship}_{ij} \le 0,\]

MP linearizes it with a big-M constant: \(\sum_j \text{Ship}_{ij} \le M \cdot \text{Open}_i\). MP takes \(M\) from the declared bounds of the variables in the expression, propagated through the expression. It does not look for bounds implied by other constraints of the model. So the same implication can be linearized with a loose or a tight constant depending only on what the declarations say, and the constant decides how strong the LP relaxation is.

The remedy is to name the expression as a variable and declare its real bound. This notebook shows the effect on a capacitated facility location model, measures it with the root (LP relaxation) bound, and looks at the constants in the model the solver actually receives.

Data: facilities with binding capacities#

import time
import numpy as np
import pandas as pd


def make_data(n_fac, n_cust, seed=3):
    rng = np.random.default_rng(seed)
    cap = rng.integers(80, 160, n_fac)
    dem = rng.integers(5, 25, n_cust)
    fix = rng.integers(300, 700, n_fac)
    xf, yf, xc, yc = (
        rng.random(n_fac),
        rng.random(n_fac),
        rng.random(n_cust),
        rng.random(n_cust),
    )
    cost = {
        (i, j): round(100 * float(np.hypot(xf[i] - xc[j], yf[i] - yc[j])), 1)
        for i in range(n_fac)
        for j in range(n_cust)
    }
    return dict(cap=cap.tolist(), dem=dem.tolist(), fix=fix.tolist(), cost=cost)


data = make_data(20, 80)
print(
    "total demand",
    sum(data["dem"]),
    " total capacity",
    sum(data["cap"]),
    " largest capacity",
    max(data["cap"]),
)

total demand 1203  total capacity 2209  largest capacity 149

Three formulations of the same model#

All three have the same integer solutions. They differ only in where the capacity information is stated.

  • implication_on_sum: the implication is written on the sum of shipments. Each Ship[i,j] is bounded by the customer’s demand, so the bound MP can derive for the sum is the total demand, far above any single capacity. The capacity is a separate constraint.

  • named_and_bounded: the sum gets a name, Out[i], whose declared upper bound is the capacity. The implication is written on Out[i].

  • named_capacity_row: Out[i] without a declared upper bound, and the capacity written as a constraint Out[i] <= cap[i].

BASE = r"""
set I; set J;
param cap {I}; param dem {J}; param fix {I}; param cost {I, J};
var Open {I} binary;
var Ship {i in I, j in J} >= 0, <= dem[j];
minimize Total: sum {i in I} fix[i] * Open[i] + sum {i in I, j in J} cost[i,j] * Ship[i,j];
subject to Demand {j in J}: sum {i in I} Ship[i,j] = dem[j];
"""

MODELS = {
    "implication_on_sum": BASE
    + r"""
subject to Capacity {i in I}: sum {j in J} Ship[i,j] <= cap[i];
subject to Use {i in I}: Open[i] = 0 ==> sum {j in J} Ship[i,j] <= 0;
""",
    "named_and_bounded": BASE
    + r"""
var Out {i in I} >= 0, <= cap[i];
subject to Out_Def {i in I}: Out[i] = sum {j in J} Ship[i,j];
subject to Use {i in I}: Open[i] = 0 ==> Out[i] <= 0;
""",
    "named_capacity_row": BASE
    + r"""
var Out {i in I} >= 0;
subject to Out_Def {i in I}: Out[i] = sum {j in J} Ship[i,j];
subject to Capacity {i in I}: Out[i] <= cap[i];
subject to Use {i in I}: Open[i] = 0 ==> Out[i] <= 0;
""",
}


def load(model, data, ampl_options=None):
    ampl = AMPL()
    for name, value in (ampl_options or {}).items():
        ampl.option[name] = value
    ampl.eval(model)
    n_fac, n_cust = len(data["cap"]), len(data["dem"])
    ampl.set["I"] = range(n_fac)
    ampl.set["J"] = range(n_cust)
    for name in ("cap", "fix"):
        ampl.param[name] = dict(enumerate(data[name]))
    ampl.param["dem"] = dict(enumerate(data["dem"]))
    ampl.param["cost"] = data["cost"]
    return ampl


def solver_objective(ampl):
    """The objective value the solver reports (last 'objective <value>' in solve_message).

    AMPL's own value, ampl.obj["Total"].value(), is recomputed from the returned
    variables. For a relaxation of a model with logical expressions (if-then, ==>)
    the two differ: AMPL evaluates the logic exactly at a fractional point."""
    import re

    found = re.findall(r"objective ([-+0-9.eE]+)", ampl.get_value("solve_message"))
    return float(found[-1]) if found else ampl.obj["Total"].value()


def solve(model, data, relax=False, solver_options="", ampl_options=None):
    """Solve once; relax=True gives the root LP bound (MP's auxiliary binaries relaxed too)."""
    ampl = load(model, data, ampl_options)
    opts = "outlev=0 " + solver_options + (" relax=1 chk:mode=0" if relax else "")
    start = time.time()
    ampl.solve(solver=SOLVER_MILP, verbose=False, **{SOLVER_MILP + "_options": opts})
    seconds = time.time() - start
    message = ampl.get_value("solve_message")
    nodes = next(
        (line.split()[0] for line in message.splitlines() if "node" in line), ""
    )
    return dict(
        status=ampl.solve_result,
        objective=solver_objective(ampl),
        nodes=nodes,
        seconds=round(seconds, 2),
    )

Compare the root bound#

The root bound is the optimum with integrality relaxed (relax=1, which also relaxes the binaries MP adds). The closer it is to the integer optimum, the less the branch-and-bound has to do.

rows = []
for name, model in MODELS.items():
    root = solve(model, data, relax=True)
    mip = solve(model, data, solver_options="timelim=120")
    rows.append(
        dict(
            formulation=name,
            root_bound=round(root["objective"], 1),
            optimum=round(mip["objective"], 1),
            status=mip["status"],
            nodes=mip["nodes"],
            seconds=mip["seconds"],
        )
    )
table = pd.DataFrame(rows).set_index("formulation")
table["root_gap_%"] = (100 * (table.optimum - table.root_bound) / table.optimum).round(
    1
)
table
root_bound optimum status nodes seconds root_gap_%
formulation
implication_on_sum 14587.8 21626.7 solved 7 0.30 32.5
named_and_bounded 19743.6 21626.7 solved 1 0.12 8.7
named_capacity_row 19743.6 21626.7 solved 1 0.12 8.7

All three reach the same optimum. The implication on the sum has a much weaker root bound, and branch-and-bound needs more nodes to close it. How much time this costs depends on the solver and the instance: a strong solver can recover part of the gap with presolve and cuts. The root gap is the solver-independent part of the story.

named_capacity_row matches named_and_bounded although its capacity is a constraint: a constraint on a single variable is a bound in disguise. AMPL’s presolve turns it into a bound, and so does MP when AMPL presolve is off:

root_row_nopresolve = solve(
    MODELS["named_capacity_row"], data, relax=True, ampl_options={"presolve": 0}
)
print(
    "named_capacity_row, AMPL presolve off: root bound",
    round(root_row_nopresolve["objective"], 1),
    " vs",
    table.loc["named_capacity_row", "root_bound"],
    "with presolve",
)

named_capacity_row, AMPL presolve off: root bound 19743.6  vs 19743.6 with presolve

implication_on_sum has a capacity constraint too, yet stays weak: there the constraint is on an expression, the sum of shipments, and a constraint on an expression never reaches the big-M sizing. So: a bound on a variable, declared or in a single-variable constraint, reaches the big-M sizing; a constraint on an expression does not. A limit on a sum, a difference or any other expression must be given its own variable with a bound.

Look at what the solver receives#

Export the solver’s model with writeprob on a tiny instance where the capacity is below the total demand. In the LP file, look for the rows with the binary Open variables (here c# columns, since the HiGHS LP writer does not keep AMPL names) and read the coefficient of the binary: that is \(-M\).

import os
import tempfile

tiny = dict(
    cap=[15, 12],
    dem=[8, 9, 8],
    fix=[620, 647],
    cost={
        (0, 0): 64.2,
        (0, 1): 5.2,
        (0, 2): 31.6,
        (1, 0): 46.7,
        (1, 1): 41.9,
        (1, 2): 42.9,
    },
)
tmpdir = tempfile.mkdtemp()
for name in ("implication_on_sum", "named_and_bounded"):
    path = os.path.join(tmpdir, name + ".lp")
    solve(MODELS[name], tiny, solver_options=f"writeprob={path}")
    lines = open(path).read().splitlines()
    rows = [l for l in lines[lines.index("st") + 1 : lines.index("bounds")]]
    print(f"--- {name}: total demand {sum(tiny['dem'])}, capacities {tiny['cap']}")
    print("\n".join(r for r in rows if "<=" in r))

--- implication_on_sum: total demand 25, capacities [15, 12]
 r3: +1 c0 +1 c1 +1 c2 <= +15
 r4: +1 c3 +1 c4 +1 c5 <= +12
 r5: +1 c0 +1 c1 +1 c2 -25 c6 <= +0
 r6: -1 c0 -1 c1 -1 c2 <= -0
 r7: +1 c3 +1 c4 +1 c5 -25 c7 <= +0
 r8: -1 c3 -1 c4 -1 c5 <= -0
--- named_and_bounded: total demand 25, capacities [15, 12]
 r5: +1 c0 -15 c2 <= +0
 r6: +1 c1 -12 c3 <= +0

With the implication on the sum, the constant multiplying each binary is the total demand, 25, the only bound MP can derive for the sum. With the named variable it is each facility’s capacity. The same check works on any model: export the solver’s model, find the rows MP generated for a logical constraint, and compare their constants with what the data says is possible.

And with no bound at all#

If a variable in a logical constraint has no finite bound, MP cannot linearize the constraint and the solve ends with a message, not an exception. Always check solve_result:

free = MODELS["implication_on_sum"].replace(
    "var Ship {i in I, j in J} >= 0, <= dem[j];", "var Ship {i in I, j in J} >= 0;"
)
ampl = load(free, data)
ampl.solve(solver=SOLVER_MILP, verbose=False)
print("solve_result:", ampl.solve_result)
print(ampl.get_value("solve_message"))

solve_result: failure
HiGHS 1.15.1: Set bounds on variables
participating in logical expressions,
or use option cvt:bigM (with caution).
See more: mp.ampl.com/modeling-numerics.html

The message also mentions cvt:bigM, a default constant for unbounded variables. Treat it as a last resort: a guessed constant is either too loose (weak relaxation) or too tight (it cuts off solutions). A declared bound that follows from the data is always better.

Takeaways#

  • MP sizes big-M constants from declared variable bounds, propagated through the expression in the logical constraint.

  • To tighten a constant, name the expression as a variable and declare the bound that follows from the data (capacity, time window, physical limit).

  • A constraint on a single variable is turned into a bound (by AMPL presolve, or by MP); a constraint on an expression stays a row and does not reach the big-M sizing.

  • Check with two cheap measurements before benchmarking: the root bound (relax=1, or option relax_integrality 1 for the model’s own integer variables), and the exported solver model (writeprob).

  • An indicator constraint, where a solver supports it natively, does not need the constant, but it is not automatically faster than a tight big-M row; measure both.

More on choosing between formulations: the MP modeling guide at https://mp.ampl.com/modeling-formulations.html, and the related notebooks under https://colab.ampl.com/tags/formulation-comparison.html.