MIP: invalid cut from a redundant variable bound gives a suboptimal solution reported as Optimal at gap 0

[EamonHetherton]

Summary

HighsTransformedLp::transform() ignores the redundant flag returned by
HighsImplications::cleanupVub() / cleanupVlb(). A redundant variable bound is deliberately
returned untightened for the caller to discard; transform() substitutes it anyway. Cut
generation is then given a false upper bound on the transformed variable, and the cut that comes
out is not globally valid. It can remove the optimum, after which the search closes and reports the
incumbent as Optimal at a zero gap, with no warning and no violated row.

Two small models that current master/latest solves incorrectly are attached as 3170-1.mps.txt (39
rows, 30 columns, 137 nonzeros, 10 binaries) and 3170-2.mps.txt (40 rows, 24 columns, 125
nonzeros, 8 binaries) — rename them to .mps.

Reproducing

3170-1.mps.txt fails with a single option and nothing else:

$ highs --presolve off 3170-1.mps

Solving report
  Model             3170-1
  Status            Optimal
  Primal bound      -21262719.4302
  Dual bound        -21262719.4302
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    -21262719.4302 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Nodes             1

The optimum is -21446579.364727, so this is suboptimal by 183,860 — a relative error of
0.86%, about 86x the default relative gap tolerance. It is not a tolerance artefact.

3170-2.mps.txt needs the option set below (with default options it is solved correctly), and then
reports -516096644.644 against an optimum of -516307714.782682 — suboptimal by 211,070, a
relative error of 0.041%, again at Gap 0% (tolerance: 0%), 1 node, zero violations:

threads=1, random_seed=0, parallel=off, run_crossover=on,
mip_rel_gap=0, mip_abs_gap=1e-6,
mip_feasibility_tolerance=1e-7, primal_feasibility_tolerance=1e-7, dual_feasibility_tolerance=1e-7,
all six mip_heuristic_run_* = false,
presolve=off

Both solve in a few hundredths of a second at 1 node. The failure is silent: the returned point is
feasible, integral and self-consistent, and the gap is certified 0.

Both optima were verified independently of the MIP solver, by enumerating all 2^n binary
assignments and solving the resulting LP for each — so the claim above does not rest on trusting
the branch-and-bound search:

import itertools, highspy

def brute_force(path):
    h = highspy.Highs(); h.setOptionValue("output_flag", False)
    h.readModel(path)
    lp = h.getLp()
    ints = [j for j in range(lp.num_col_)
            if lp.integrality_[j] == highspy.HighsVarType.kInteger]
    best = float("inf")
    for bits in itertools.product((0.0, 1.0), repeat=len(ints)):
        g = highspy.Highs(); g.setOptionValue("output_flag", False)
        g.readModel(path)
        for j, v in zip(ints, bits):
            g.changeColBounds(j, v, v)
            g.changeColIntegrality(j, highspy.HighsVarType.kContinuous)
        g.run()
        if g.getModelStatus() == highspy.HighsModelStatus.kOptimal:
            best = min(best, g.getObjectiveValue())
    return best

print(brute_force("3170-1.mps"))   # -21446579.364727
print(brute_force("3170-2.mps"))   # -516307714.782682

On the presolve option. Disabling presolve is only what makes these two models expose the
defect — presolve itself is not doing anything wrong, and the faulty code is in cut generation,
which runs either way. Disabling it changes the relaxation that reaches the separators, which is
what determines whether the unsound substitution below is reached at all.

Cause

highs/mip/HighsTransformedLp.cpp, transform().

Cut generation substitutes bounds to move each variable into [0, ub-lb]. When it substitutes a
variable upper bound x <= a*z + b (z binary), the transformed variable is the slack
y = a*z + b - x, and the code passes upper[j] = ub - lb to the cut generator as y's upper
bound. That is sound only if the variable bound is at least as tight as the simple bound, i.e.
maxValue() <= ub. The code checks exactly that and tries to repair it:

if (bestVub[col].first != -1 &&
    bestVub[col].second.maxValue() > ub + mip.mipdata_->feastol) {
  bool redundant = false;
  bool infeasible = false;
  mip.mipdata_->implications.cleanupVub(col, bestVub[col].first,
                                        bestVub[col].second, ub, redundant,
                                        infeasible, false);
}

but the redundant output is discarded. cleanupVub() tightens the bound only when it is not
redundant:

if (minub >= ub - mipsolver.mipdata_->feastol) {
  redundant = true;          // left unchanged, for the caller to discard
} else if (maxub > ub + mipsolver.mipdata_->epsilon) {
  ...                        // tightened so that maxValue() == ub
}

so on the redundant path the variable bound comes back with maxValue() still far above ub,
and the invariant transform() just tried to establish still fails. cleanupVarbounds() handles
this correctly — it collects redundant bounds and erases them — but transform() simply carries
on and substitutes it. cleanupVlb() has the identical contract and the identical caller bug.

On a larger model where I first hit this, a continuous column with global bounds
[0.0298, 0.0683] (so ub-lb = 0.0385) was substituted with

x <= 0.29538846115832368 * z + 0.4420490388416824       minValue 0.4421   maxValue 0.7374

cleanupVub() returned redundant = 1 and changed nothing. The real range of the slack is
[0, 0.7077], but the cut generator was told y <= 0.0385 — 18x too small. The cut that came out
had a minimum activity over the global box of +1.28 against a right-hand side of -0.55: it is
violated by every point in the box. Once such a cut is in the global pool every node LP is
infeasible, the tree closes at once, and whatever incumbent existed is declared optimal.

That last detail may also be worth checking against reports where a feasible MIP is returned
Infeasible — a globally infeasible cut added before any incumbent exists produces exactly that.

How the state arises

addVUB() refuses to store an already-redundant bound, so this needs one that was valid when
stored and became redundant later:

1. probing stores a variable upper bound for x with minValue = c and maxValue = ub(x)
   (probing-derived bounds always have maxValue equal to the column's bound at derivation time);
2. x's upper bound is later tightened to <= c — for instance because the binary gets fixed;
3. getBestVub() had already selected that bound for the current separation round, and
   transform() then takes the redundant path above.

Step 3 is why the binary being fixed does not save you: getBestVub() skips fixed binaries, but
the selection happened earlier in the round.

Fix

A candidate fix, a regression test using the two attached models that fails without it, and the
supporting measurements are on PR #3172, so as to keep this report to the diagnosis.

Version

master @ 04024d7 and latest @ bc30ccf (1.15.1). Also reproduces on 1.13.1; I have not tried other versions systematically.

Related

There is a second, independent unsoundness in the same area — cut coefficients built with one bound
substitution and untransformed with another: #3171.

Both come from cut construction reading bound-substitution state back out of a shared, still-
mutating HighsTransformedLp, rather than the cut carrying the transformation it was built with.
They are independent in practice though: neither fix subsumes the other, and each can be reviewed
on its own.

#3173 is a third defect found on the same class of model, but in presolve rather than in cut
generation - it is decided before branch-and-bound runs and is unaffected by the fix here.

[jajhall]

See #3172 @Opt-Mucca and @fwesselm

[Opt-Mucca]

This is a bug. Good find, and a nice description. Normally the larger bound would simply lead to weaker cuts, but everything gets confusing after complementing and creating covers etc. I'll leave another comment in the PR.