Fix non-rigorous crossed bound results

pyomo/contrib/gdpopt/algorithm_base_class.py

@@ -14,6 +14,7 @@
 from pyomo.common.errors import DeveloperError
 from pyomo.common.modeling import unique_component_name
 from pyomo.contrib.gdpopt.config_options import _add_common_configs
+from pyomo.contrib.gdpopt.convexity import model_is_not_certified_convex
 from pyomo.contrib.gdpopt.create_oa_subproblems import (
     add_util_block,
     add_disjunct_list,
@@ -37,6 +38,7 @@
 class _GDPoptAlgorithm:
     CONFIG = ConfigBlock("GDPopt")
     _add_common_configs(CONFIG)
+    _requires_model_convexity = False
 
     def __init__(self, **kwds):
         """
@@ -57,6 +59,8 @@ def __init__(self, **kwds):
 
         self.incumbent_boolean_soln = None
         self.incumbent_continuous_soln = None
+        self._bounds_crossed_without_certified_convexity = False
+        self._model_is_not_certified_convex = False
 
         self.original_obj = None
         self._dummy_obj = None
@@ -113,6 +117,8 @@ def solve(self, model, **kwds):
 
         config = self.config(kwds.pop('options', {}), preserve_implicit=True)
         config.set_value(kwds)
+        self._bounds_crossed_without_certified_convexity = False
+        self._model_is_not_certified_convex = False
 
         with lower_logger_level_to(config.logger, tee=config.tee):
             self._log_solver_intro_message(config)
@@ -207,6 +213,10 @@ def _gather_problem_info_and_solve_non_gdps(self, model, config):
         logger = config.logger
 
         self._create_pyomo_results_object_with_problem_info(model, config)
+        self._model_is_not_certified_convex = (
+            self._requires_model_convexity
+            and model_is_not_certified_convex(model, config.eigenvalue_tolerance)
+        )
         # Check if this problem actually has any discrete decisions. If not,
         # just solve it.
         problem = self.pyomo_results.problem
@@ -372,6 +382,17 @@ def bounds_converged(self, config):
                 self._load_infeasible_termination_status(config)
             elif self.LB == float('-inf') and self.UB == float('-inf'):
                 self._load_infeasible_termination_status(config)
+            elif (
+                self.LB > self.UB
+                and self._requires_model_convexity
+                and self._model_is_not_certified_convex
+            ):
+                self._bounds_crossed_without_certified_convexity = True
+                self._log_current_state(config.logger, '')
+                config.logger.info(
+                    'GDPopt exiting--bounds crossed without certified model convexity.'
+                )
+                self.pyomo_results.solver.termination_condition = tc.feasible
             else:
                 # if they've crossed, then the gap is actually 0: Update the
                 # dual (discrete problem) bound to be equal to the primal
@@ -518,8 +539,16 @@ def _get_final_pyomo_results_object(self):
         """
         results = self.pyomo_results
         # Finalize results object
-        results.problem.lower_bound = self.LB
-        results.problem.upper_bound = self.UB
+        if self._bounds_crossed_without_certified_convexity:
+            if self.objective_sense is minimize:
+                results.problem.lower_bound = float('-inf')
+                results.problem.upper_bound = self.UB
+            else:
+                results.problem.lower_bound = self.LB
+                results.problem.upper_bound = float('inf')
+        else:
+            results.problem.lower_bound = self.LB
+            results.problem.upper_bound = self.UB
         results.solver.iterations = self.iteration
         results.solver.timing = self.timing
         results.solver.user_time = self.timing.total

pyomo/contrib/gdpopt/config_options.py

@@ -529,6 +529,14 @@ def _add_tolerance_configs(CONFIG):
             description="Tolerance for bound convergence.",
         ),
     )
+    CONFIG.declare(
+        "eigenvalue_tolerance",
+        ConfigValue(
+            default=1e-10,
+            domain=NonNegativeFloat,
+            description=("Numerical tolerance for eigenvalue-based convexity checks."),
+        ),
+    )
 
 
 def _add_ldsda_configs(CONFIG):

pyomo/contrib/gdpopt/convexity.py

@@ -0,0 +1,117 @@
+# ____________________________________________________________________________________
+#
+# Pyomo: Python Optimization Modeling Objects
+# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC
+# Under the terms of Contract DE-NA0003525 with National Technology and Engineering
+# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this
+# software.  This software is distributed under the 3-clause BSD License.
+# ____________________________________________________________________________________
+
+"""Conservative convexity checks for outer approximation algorithms.
+
+Outer approximation algorithms that linearize the nonlinear constraints at trial
+points (GDPopt LOA, MindtPy OA, MindtPy ECP) only produce a valid relaxation of
+the original problem when that problem is convex. Applied to a nonconvex problem,
+the linearizations can cut off feasible points, so the resulting "dual bound" is
+not a rigorous bound and must not be reported as certifying global optimality.
+
+Algorithms that build their relaxation from McCormick envelopes instead (GDPopt
+GLOA, MindtPy GOA) do produce a valid relaxation for nonconvex problems, so their
+bounds are rigorous and are unaffected by this module.
+
+The detection here is deliberately conservative: anything this module cannot
+positively certify as convex is reported as not certified convex.
+"""
+
+from pyomo.common.dependencies import numpy as np
+from pyomo.core import Block, Constraint, Objective, minimize, value
+from pyomo.core.base.enums import SortComponents
+from pyomo.gdp import Disjunct
+from pyomo.repn.quadratic import QuadraticRepnVisitor
+from pyomo.repn.util import OrderedVarRecorder
+
+
+def _quadratic_matrix(repn):
+    """Build the symmetric matrix Q of a quadratic repn.
+
+    Returns None if any quadratic coefficient cannot be evaluated to a number.
+    """
+    var_to_idx = {}
+    for var_ids in repn.quadratic:
+        for var_id in var_ids:
+            if var_id not in var_to_idx:
+                var_to_idx[var_id] = len(var_to_idx)
+
+    q_matrix = np.zeros((len(var_to_idx), len(var_to_idx)))
+    for (var_id1, var_id2), coef in repn.quadratic.items():
+        coef_val = value(coef, exception=False)
+        if coef_val is None:
+            return None
+        idx1 = var_to_idx[var_id1]
+        idx2 = var_to_idx[var_id2]
+        if var_id1 == var_id2:
+            q_matrix[idx1][idx1] += coef_val
+        else:
+            half_coef = 0.5 * coef_val
+            q_matrix[idx1][idx2] += half_coef
+            q_matrix[idx2][idx1] += half_coef
+
+    return q_matrix
+
+
+def quadratic_curvature(expr, eigenvalue_tolerance):
+    """Classify the curvature of a quadratic expression.
+
+    Returns 1 if the quadratic form is positive semidefinite (convex), -1 if it
+    is negative semidefinite (concave), 0 if it has no quadratic terms or the
+    quadratic form vanishes, and None if the curvature could not be determined.
+    """
+    recorder = OrderedVarRecorder({}, {}, SortComponents.deterministic)
+    repn = QuadraticRepnVisitor({}, var_recorder=recorder).walk_expression(expr)
+    if repn.nonlinear is not None:
+        return None
+    if repn.quadratic is None:
+        return 0
+
+    q_matrix = _quadratic_matrix(repn)
+    if q_matrix is None:
+        return None
+
+    eigenvalues = np.linalg.eigvalsh(q_matrix)
+    is_psd = all(eigenvalue >= -eigenvalue_tolerance for eigenvalue in eigenvalues)
+    is_nsd = all(eigenvalue <= eigenvalue_tolerance for eigenvalue in eigenvalues)
+    if is_psd and is_nsd:
+        return 0
+    if is_psd:
+        return 1
+    if is_nsd:
+        return -1
+    return None
+
+
+def model_is_not_certified_convex(model, eigenvalue_tolerance):
+    """Return True unless this model can be certified as convex.
+
+    A True result means an outer approximation dual bound computed for the model
+    must not be treated as rigorous.
+    """
+    for obj in model.component_data_objects(Objective, active=True, descend_into=True):
+        curvature = quadratic_curvature(obj.expr, eigenvalue_tolerance)
+        if obj.sense is minimize and curvature not in (0, 1):
+            return True
+        elif obj.sense is not minimize and curvature not in (0, -1):
+            return True
+
+    for constr in model.component_data_objects(
+        Constraint, active=True, descend_into=(Block, Disjunct)
+    ):
+        curvature = quadratic_curvature(constr.body, eigenvalue_tolerance)
+        if curvature == 0:
+            continue
+        if constr.equality:
+            return True
+        if constr.has_ub() and curvature not in (0, 1):
+            return True
+        if constr.has_lb() and curvature not in (0, -1):
+            return True
+    return False

pyomo/contrib/gdpopt/loa.py

@@ -74,6 +74,7 @@ class GDP_LOA_Solver(_GDPoptAlgorithm, _OAAlgorithmMixIn):
     _add_tolerance_configs(CONFIG)
 
     algorithm = 'LOA'
+    _requires_model_convexity = True
 
     # Override solve() to customize the docstring for this solver
     @document_kwargs_from_configdict(CONFIG, doc=_GDPoptAlgorithm.solve.__doc__)

pyomo/contrib/gdpopt/tests/test_gdpopt.py

@@ -12,14 +12,19 @@
 from contextlib import redirect_stdout
 from io import StringIO
 import logging
+from pyomo.common.timing import default_timer
 from math import fabs
 from os.path import join, normpath
 
 import pyomo.common.unittest as unittest
 from pyomo.common.log import LoggingIntercept
 from pyomo.common.collections import Bunch
 from pyomo.common.config import ConfigDict, ConfigValue
+from pyomo.common.dependencies import numpy_available
 from pyomo.common.fileutils import import_file, PYOMO_ROOT_DIR
+from pyomo.contrib.gdpopt.convexity import model_is_not_certified_convex
+from pyomo.contrib.gdpopt.gloa import GDP_GLOA_Solver
+from pyomo.contrib.gdpopt.loa import GDP_LOA_Solver
 from pyomo.contrib.gdpopt.create_oa_subproblems import (
     add_util_block,
     add_disjunct_list,
@@ -44,13 +49,14 @@
     RangeSet,
     TransformationFactory,
     SolverFactory,
+    sin,
     sqrt,
     value,
     Var,
 )
 from pyomo.gdp import Disjunct, Disjunction
 from pyomo.gdp.tests import models
-from pyomo.opt import TerminationCondition
+from pyomo.opt import SolverResults, TerminationCondition
 
 exdir = normpath(join(PYOMO_ROOT_DIR, 'examples', 'gdp'))
 
@@ -78,6 +84,123 @@
 class TestGDPoptUnit(unittest.TestCase):
     """Real unit tests for GDPopt"""
 
+    def _setup_crossed_bound_state(self, solver, sense=maximize):
+        solver.pyomo_results = SolverResults()
+        solver.pyomo_results.problem.sense = sense
+        solver.LB = 1011.6577899409375
+        solver.UB = 1000.0
+        solver.iteration = 1
+        solver._model_is_not_certified_convex = True
+        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(__name__)
+        return config
+
+    def test_loa_crossed_bounds_are_not_reported_as_global_optimal(self):
+        solver = GDP_LOA_Solver()
+        config = self._setup_crossed_bound_state(solver)
+
+        self.assertTrue(solver.bounds_converged(config))
+        results = solver._get_final_pyomo_results_object()
+
+        self.assertIs(
+            results.solver.termination_condition, TerminationCondition.feasible
+        )
+        self.assertEqual(results.problem.lower_bound, 1011.6577899409375)
+        self.assertEqual(results.problem.upper_bound, float('inf'))
+
+    def test_loa_nonpolynomial_model_is_not_certified_convex(self):
+        m = ConcreteModel()
+        m.x = Var(bounds=(0, 4), initialize=0.2)
+        m.y = Var(bounds=(-2, 2), initialize=0)
+        m.nonconvex = Constraint(expr=m.y <= sin(3 * m.x) + 0.2 * m.x)
+        m.choose_region = Disjunction(expr=[[m.x <= 2], [m.x >= 2]], xor=True)
+        m.obj = Objective(expr=m.y, sense=maximize)
+
+        self.assertTrue(model_is_not_certified_convex(m, 1e-10))
+
+    @unittest.skipUnless(numpy_available, 'NumPy is not available')
+    def test_loa_distinguishes_basic_quadratic_convexity(self):
+        convex = ConcreteModel()
+        convex.x = Var(bounds=(-2, 2))
+        convex.y = Var(bounds=(-2, 2))
+        convex.c = Constraint(expr=convex.x**2 + convex.y**2 <= 1)
+        convex.obj = Objective(expr=convex.x**2 + convex.y)
+
+        nonconvex = ConcreteModel()
+        nonconvex.x = Var(bounds=(-2, 2))
+        nonconvex.y = Var(bounds=(-2, 2))
+        nonconvex.c = Constraint(expr=nonconvex.x * nonconvex.y <= 1)
+        nonconvex.obj = Objective(expr=nonconvex.y)
+
+        self.assertFalse(model_is_not_certified_convex(convex, 1e-10))
+        self.assertTrue(model_is_not_certified_convex(nonconvex, 1e-10))
+
+    @unittest.skipUnless(numpy_available, 'NumPy is not available')
+    def test_loa_certifies_psd_quadratic_cross_terms(self):
+        m = ConcreteModel()
+        m.x = Var(bounds=(-2, 2))
+        m.y = Var(bounds=(-2, 2))
+        m.c = Constraint(expr=m.x**2 + m.x * m.y + m.y**2 <= 4)
+        m.obj = Objective(expr=m.x)
+
+        self.assertFalse(model_is_not_certified_convex(m, 1e-10))
+
+    @unittest.skipUnless(numpy_available, 'NumPy is not available')
+    def test_loa_eigenvalue_tolerance_is_configurable(self):
+        m = ConcreteModel()
+        m.x = Var()
+        m.c = Constraint(expr=-1e-11 * m.x**2 <= 1)
+        m.obj = Objective(expr=m.x)
+
+        config = GDP_LOA_Solver.CONFIG()
+        self.assertFalse(model_is_not_certified_convex(m, config.eigenvalue_tolerance))
+        config.eigenvalue_tolerance = 1e-12
+        self.assertTrue(model_is_not_certified_convex(m, config.eigenvalue_tolerance))
+
+    def test_gloa_crossed_bounds_preserve_certified_optimal_behavior(self):
+        solver = GDP_GLOA_Solver()
+        config = self._setup_crossed_bound_state(solver)
+
+        self.assertTrue(solver.bounds_converged(config))
+        results = solver._get_final_pyomo_results_object()
+
+        self.assertIs(
+            results.solver.termination_condition, TerminationCondition.optimal
+        )
+        self.assertEqual(results.problem.lower_bound, 1011.6577899409375)
+        self.assertEqual(results.problem.upper_bound, 1011.6577899409375)
+
+    @unittest.skipUnless(numpy_available, 'NumPy is not available')
+    @unittest.skipUnless(LOA_solvers_available, 'Required subsolvers are not available')
+    def test_loa_nonconvex_solve_reports_crossed_bounds_as_feasible(self):
+        m = ConcreteModel()
+        m.x = Var(bounds=(0, 4), initialize=0.2)
+        m.y = Var(bounds=(-2, 2), initialize=0)
+        m.nonconvex = Constraint(expr=m.y <= sin(3 * m.x) + 0.2 * m.x)
+        m.left = Disjunct()
+        m.right = Disjunct()
+        m.left.region = Constraint(expr=m.x <= 2)
+        m.right.region = Constraint(expr=m.x >= 2)
+        m.choose_region = Disjunction(expr=[m.left, m.right], xor=True)
+        m.obj = Objective(expr=m.y, sense=maximize)
+
+        results = SolverFactory('gdpopt.loa').solve(
+            m,
+            mip_solver=mip_solver,
+            nlp_solver=nlp_solver,
+            init_algorithm='set_covering',
+        )
+
+        self.assertIs(
+            results.solver.termination_condition, TerminationCondition.feasible
+        )
+        self.assertEqual(results.problem.lower_bound, value(m.obj))
+        self.assertEqual(results.problem.upper_bound, float('inf'))
+
     @unittest.skipUnless(
         SolverFactory(mip_solver).available(), "MIP solver not available"
     )