Adding compatibility checks for HiGHS (and perhaps other solvers?)

[tvalenciaz]

Writing and running some examples with HiGHS, I identified a couple of issues where we could check compatibility from the start. I'm not sure if there was already a discussion about whether it is our job or the user's to verify that the solver they intend to use does support the features requested, but we could issue more helpful errors:

1. If using HiGHS with PH and problem has integers/binaries. HiGHS does not support MIQP and just fails with status= unknown, TerminationCondition=unknown. It's true that if the user runs with --tee, then the solver does print "cannot solve MIQP with HiGHS". But we could make it easier for the user requiring --linearize-proximal-terms , or at least issue a more descriptive error.
2. HiGHS does not return reduced costs*, so --reduced-cost-spoke and related features cannot be used.
   *I am not 100% sure about this, but I have noticed that the suffix attribute that is supposed to hold the reduced costs is empty when I use HiGHS and the error thrown is a bit cryptic. My hypothesis is that the solver does not return them, so we could issue a more helpful error saying that. This was tested with the hydro example in generic_tester.py as well as with the "relief" example introduced in Adding (basic) LShaped support to generic cylinders #761

[bknueven]

For (1), I think we should print a message if the solver fails at PH iteration 1 and we're using quadratic proximal terms. Something like "solver failed at PH iteration 1, likely because solver {solver} does not support quadratic objectives". Maybe if we got fancier we could try to extract the error message and re-raise it -- not sure where that message is coming from in the case of HiGHS.

For (2), this sounds like something that should be investigated and fixed in Pyomo. HiGHS surely reports reduced costs, Pyomo just isn't extracting them on request. Or, this could be an issue with mpi-sppy and the new Pyomo solver interfaces. I would suggest we investigate the root cause before putting in a HiGHS-specific patch.

[DLWoodruff]

Investigated part 2 (reduced costs with HiGHS). The root cause is not HiGHS — HiGHS does report reduced costs. It's a bug in Pyomo's APPSI legacy solver wrapper that surfaces because of how mpi-sppy loads solutions, and it affects every APPSI solver (appsi_gurobi fails identically).

Reproduction (LP with an rc IMPORT suffix)

Path	primals	duals	reduced costs
appsi_highs.solve(m) (default load_solutions=True)	ok	ok	ok (rc[x]=-1.0)
appsi_highs.solve(m, load_solutions=False) + m.solutions.load_from(results)	ok	ok	empty (len(rc)=0)
same with gurobi/cplex shell solvers	ok	ok	ok

mpi-sppy always uses the second path (spopt.py solves with load_solutions=False, then s.solutions.load_from(results), so it can inspect the termination condition before loading). That path silently drops reduced costs for any APPSI solver — confirming this is "an issue with mpi-sppy and the new Pyomo solver interfaces," not a HiGHS limitation.

Root cause — a Pyomo bug

pyomo/contrib/appsi/base.py, LegacySolverInterface.solve, the load_solutions=False branch:

# duals — CORRECT (per-constraint, keyed by the constraint symbol)
legacy_soln.constraint[symbol_map.getSymbol(c)] = {'Dual': val}
...
# reduced costs — BUG
if hasattr(model, 'rc') and model.rc.import_enabled():
    for v, val in results.solution_loader.get_reduced_costs().items():
        legacy_soln.variable['Rc'] = val          # wrong key, scalar, overwritten each iteration

It should attach 'Rc' onto each variable's own entry (which the primal loop already created as {'Value': ...}), i.e. legacy_soln.variable[symbol_map.getSymbol(v)]['Rc'] = val — mirroring the dual handling. PyomoModel.select() lowercases 'Rc'→'rc' and routes per-variable entries to the rc suffix, but the bogus literal 'Rc' key never maps to a variable, so the suffix stays empty. That is exactly why duals load but reduced costs don't. This looks worth reporting upstream to Pyomo.

What the user currently sees

With --reduced-cost-spoke, the empty rc suffix trips _assert_consensus_rc_loaded in reduced_costs_spoke.py, raising RuntimeError: reduced cost not loaded for {v.name}; vars_to_load did not cover the consensus group — which blames vars_to_load and is misleading; the real cause is the dropped suffix.

A workaround is available (if we want one in mpi-sppy)

Even after load_solutions=False, the LegacySolver wrapper retains the solution internally, and s._solver_plugin.get_reduced_costs(vars_to_load=...) returns correct values. The spoke already special-cases old-style persistent solvers (load_rc); we could add a parallel branch for the new APPSI solvers (which sputils.is_persistent does not currently recognize, since it checks only the old PersistentSolver base class).

Environment: pyomo 6.9.1, highspy 1.14.0.

[DLWoodruff]

I'm going to try to submit a Pyomo PR to fix this.

[DLWoodruff]

Filed the upstream fix as Pyomo/pyomo#3978 — it corrects the reduced-cost loading bug in both legacy solver wrappers (APPSI and the new contrib.solver interface) so rc populates correctly on the load_solutions=False + load_from path that mpi-sppy uses. Once that lands and is in a Pyomo release, --reduced-cost-spoke with HiGHS (and any APPSI solver) will work without a HiGHS-specific patch on our side.

[DLWoodruff]

Part 1 (quadratic proximal term / HiGHS MIQP) is addressed in #768. It adds a hybrid check in phbase.iterk_loop: a proactive has_capability('quadratic_objective') probe (catches solvers like cbc/glpk before solving) plus a reactive "no solution anywhere at iteration 1" check via allreduce_or for solvers that expose no capability info (HiGHS via the APPSI / contrib.solver interfaces). The error points the user at --linearize-proximal-terms. So both parts now have PRs: part 1 → #768, part 2 (reduced costs) → Pyomo/pyomo#3978.