GDPopt LOA reports non-rigorous nonconvex bounds as optimal LB=UB

[bernalde]

Summary

GDPopt LOA can report TerminationCondition.optimal with problem.lower_bound == problem.upper_bound on nonconvex GDPs even when those bounds are not rigorous.

This is not a request for LOA to prove global optimality on nonconvex problems. Pyomo already documents that "For nonconvex problems, LOA may not report rigorous lower/upper bounds." The issue is that the returned SolverResults object still uses the unqualified global-looking fields:

• results.solver.termination_condition = optimal
• results.problem.lower_bound = <local incumbent>
• results.problem.upper_bound = <same local incumbent>

For benchmark and downstream tooling, that is indistinguishable from a globally certified solve unless the caller already knows the model is nonconvex and knows this GDPopt LOA caveat.

Minimal result-finalization reproducer

This reproducer does not require external solvers. It isolates the relevant result semantics: when a LOA run has a local primal incumbent and a non-rigorous OA/discrete bound that has crossed it, the generic GDPopt convergence logic force-updates the dual bound to the primal bound and returns optimal with equal bounds.

import logging

import pyomo.environ as pyo
from pyomo.common.timing import default_timer
from pyomo.contrib.gdpopt.loa import GDP_LOA_Solver
from pyomo.opt import SolverResults

solver = GDP_LOA_Solver()
solver.pyomo_results = SolverResults()
solver.pyomo_results.problem.sense = pyo.maximize

# This mimics a nonconvex LOA state:
# - LB is a feasible incumbent from a local NLP subproblem.
# - UB is the OA/discrete bound. For a nonconvex problem this bound is not
#   necessarily globally valid, and it has crossed the incumbent.
solver.LB = 1011.6577899409375
solver.UB = 1000.0
solver.iteration = 1
solver.timing.main_timer_start_time = default_timer()
solver.timing.total = 0.0

config = solver.CONFIG()
config.bound_tolerance = 1e-6
config.logger = logging.getLogger("gdpopt-loa-mwe")
config.logger.handlers[:] = [logging.StreamHandler()]
config.logger.setLevel(logging.INFO)

print("before", solver.LB, solver.UB, solver.pyomo_results.solver.termination_condition)
print("converged", solver.bounds_converged(config))
print("after", solver.LB, solver.UB, solver.pyomo_results.solver.termination_condition)

results = solver._get_final_pyomo_results_object()
print(
    "results",
    results.problem.lower_bound,
    results.problem.upper_bound,
    results.solver.termination_condition,
)

Observed output with Pyomo 6.10.0:

        1                      1011.65779    1011.65779      0.00%      0.00
GDPopt exiting--bounds have converged or crossed.
before 1011.6577899409375 1000.0 unknown
converged True
after 1011.6577899409375 1011.6577899409375 optimal
results 1011.6577899409375 1011.6577899409375 optimal

The problematic behavior is not that the internal crossed-bound logic exists for algorithms with rigorous bounds. The problematic behavior is that LOA uses the same finalization path even though its dual bound can be non-rigorous on nonconvex models.

Downstream evidence from GDPLib

This surfaced while benchmarking GDPLib's HDA model:

• GDPLib HDA issue: Refactor hda to fix benchmark issues SECQUOIA/gdplib#67
• GDPLib PR: Fix hda hull and GLOA benchmark paths SECQUOIA/gdplib#130
• GDPLib central benchmark tracker: Central tracker: model benchmark failures across GDP strategies SECQUOIA/gdplib#60

For the HDA instance, the local LOA benchmark row reported:

strategy: gdpopt.loa
solver profile: gams-local
termination: optimal
lower bound: 1011.6577899409375
upper bound: 1011.6577899409375
iterations: 1

However, the HDA model is nonconvex, and other GDPopt/global-capable routes found a better feasible solution:

gdpopt.gloa + gams-local: optimal, objective 1105.3441871686007
gdpopt.ric + gams-local: optimal, objective 1105.3441871686007
gdpopt.ric + gams-gurobi: optimal, objective 1105.3441871686007
gdpopt.ric + gams-baron: optimal, objective 1105.339918385673

For a maximization problem, a reported LOA upper bound of 1011.6577899409375 is therefore not a valid global upper bound once a feasible point around 1105.34 is known. The LB = UB result gives the appearance of a globally certified incumbent even though it is a local/non-rigorous LOA result.

Expected behavior

For LOA on models where the OA/discrete bound is not known to be rigorous, GDPopt should avoid returning globally certified-looking result metadata.

Possible approaches:

• Return a non-global termination condition such as feasible or locallyOptimal when LOA converges only by non-rigorous/crossed bounds on a nonconvex problem.
• Preserve the incumbent as the primal bound but avoid force-setting the dual bound equal to the incumbent when the dual bound is not rigorous.
• Add explicit result metadata indicating that the reported bound is not rigorous.
• Add an option such as assume_convex=True/False or certify_bounds=True/False, with conservative default result semantics for LOA unless convexity/global validity is asserted.

The important downstream requirement is that JSON/benchmark consumers should not have to infer from solver.name == "GDPopt ... - LOA" that optimal and LB = UB may not mean globally optimal.

Environment

Observed locally with:

Python 3.12.13
Pyomo 6.10.0
OS/platform: Linux WSL2, glibc 2.39

[bernalde]

I added a full model MWE here, in addition to the solver-free result-finalization reproducer in the issue body.

This example uses open solvers (glpk for the MILP master and ipopt for NLP subproblems). It is intentionally small: the reported LOA upper bound is invalidated by an explicit feasible point in the other disjunct.

import math

import pyomo.environ as pyo
from pyomo.gdp import Disjunct, Disjunction


def build_model():
    m = pyo.ConcreteModel()
    m.x = pyo.Var(bounds=(0, 4), initialize=0.2)
    m.y = pyo.Var(bounds=(-2, 2), initialize=0)

    # Nonconvex feasible region: y <= sin(3*x) + 0.2*x.
    m.nonconvex = pyo.Constraint(expr=m.y <= pyo.sin(3 * m.x) + 0.2 * m.x)

    # A small GDP structure. The better feasible point is in the right branch,
    # but LOA terminates after the left-branch local solution.
    m.left = Disjunct()
    m.right = Disjunct()
    m.left.region = pyo.Constraint(expr=m.x <= 2)
    m.right.region = pyo.Constraint(expr=m.x >= 2)
    m.choose_region = Disjunction(expr=[m.left, m.right], xor=True)

    m.obj = pyo.Objective(expr=m.y, sense=pyo.maximize)
    return m


model = build_model()
results = pyo.SolverFactory("gdpopt.loa").solve(
    model,
    mip_solver="glpk",
    nlp_solver="ipopt",
    init_algorithm="set_covering",
    bound_tolerance=1e-6,
    tee=False,
)

print("LOA termination:", results.solver.termination_condition)
print("LOA lower bound:", results.problem.lower_bound)
print("LOA upper bound:", results.problem.upper_bound)
print("LOA incumbent x:", pyo.value(model.x))
print("LOA incumbent y:", pyo.value(model.y))
print("left indicator:", model.left.indicator_var.value)
print("right indicator:", model.right.indicator_var.value)

# An explicit feasible point in the right disjunct. This is enough to show that
# the LOA-reported upper bound is not a valid global upper bound.
x_feas = 2.6404
y_feas = math.sin(3 * x_feas) + 0.2 * x_feas
print("explicit feasible x:", x_feas)
print("explicit feasible y:", y_feas)
print("feasible y exceeds reported upper bound:", y_feas > results.problem.upper_bound)

Observed output with Pyomo 6.10.0:

LOA termination: optimal
LOA lower bound: 1.1069428089835942
LOA upper bound: 1.1069428089835942
LOA incumbent x: 0.5458374919734859
LOA incumbent y: 1.1069428089835942
left indicator: True
right indicator: False
explicit feasible x: 2.6404
explicit feasible y: 1.5258216961360678
feasible y exceeds reported upper bound: True

The explicit point x = 2.6404, y = sin(3*x) + 0.2*x = 1.5258216961360678 satisfies the right disjunct (x >= 2) and the nonlinear constraint. For this maximization model, that means the reported LOA upper bound 1.1069428089835942 is not a valid global upper bound. The issue is therefore not only that LOA found a local solution, but that the returned SolverResults fields present it as optimal with LB = UB.

Related downstream GDPLib links:

• GDPLib HDA issue: Refactor hda to fix benchmark issues SECQUOIA/gdplib#67
• GDPLib HDA PR: Fix hda hull and GLOA benchmark paths SECQUOIA/gdplib#130
• GDPLib central benchmark tracker: Central tracker: model benchmark failures across GDP strategies SECQUOIA/gdplib#60
• GDPLib PR note for this LOA caveat: Fix hda hull and GLOA benchmark paths SECQUOIA/gdplib#130 (comment)