Fix non-rigorous crossed bound results

pyomo/contrib/gdpopt/algorithm_base_class.py

@@ -37,6 +37,7 @@
 class _GDPoptAlgorithm:
     CONFIG = ConfigBlock("GDPopt")
     _add_common_configs(CONFIG)
+    _crossed_bounds_are_certified = True
 
     def __init__(self, **kwds):
         """
@@ -57,6 +58,8 @@ def __init__(self, **kwds):
 
         self.incumbent_boolean_soln = None
         self.incumbent_continuous_soln = None
+        self._nonrigorous_crossed_bounds = False
+        self._nonrigorous_dual_bound_possible = False
 
         self.original_obj = None
         self._dummy_obj = None
@@ -113,6 +116,8 @@ def solve(self, model, **kwds):
 
         config = self.config(kwds.pop('options', {}), preserve_implicit=True)
         config.set_value(kwds)
+        self._nonrigorous_crossed_bounds = False
+        self._nonrigorous_dual_bound_possible = False
 
         with lower_logger_level_to(config.logger, tee=config.tee):
             self._log_solver_intro_message(config)
@@ -207,6 +212,9 @@ 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._nonrigorous_dual_bound_possible = (
+            self._problem_may_have_nonrigorous_dual_bound(model)
+        )
         # Check if this problem actually has any discrete decisions. If not,
         # just solve it.
         problem = self.pyomo_results.problem
@@ -352,6 +360,9 @@ def _update_dual_bound_to_infeasible(self):
         else:
             self._update_bounds(dual=float('-inf'))
 
+    def _problem_may_have_nonrigorous_dual_bound(self, model):
+        return False
+
     def _update_primal_bound_to_unbounded(self, config):
         if self.objective_sense == minimize:
             self._update_bounds(primal=float('-inf'))
@@ -372,6 +383,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 not self._crossed_bounds_are_certified
+                and self._nonrigorous_dual_bound_possible
+            ):
+                self._nonrigorous_crossed_bounds = True
+                self._log_current_state(config.logger, '')
+                config.logger.info(
+                    'GDPopt exiting--bounds crossed without a certified dual bound.'
+                )
+                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 +540,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._nonrigorous_crossed_bounds:
+            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/loa.py

@@ -8,10 +8,11 @@
 # ____________________________________________________________________________________
 
 from collections import namedtuple
-from math import copysign
+from math import copysign, sqrt
 
 from pyomo.common.collections import ComponentMap
 from pyomo.common.config import document_kwargs_from_configdict
+from pyomo.common.dependencies import numpy as np, numpy_available
 from pyomo.common.modeling import unique_component_name
 from pyomo.contrib.gdpopt.algorithm_base_class import _GDPoptAlgorithm
 from pyomo.contrib.gdpopt.config_options import (
@@ -74,12 +75,151 @@ class GDP_LOA_Solver(_GDPoptAlgorithm, _OAAlgorithmMixIn):
     _add_tolerance_configs(CONFIG)
 
     algorithm = 'LOA'
+    _crossed_bounds_are_certified = False
 
     # Override solve() to customize the docstring for this solver
     @document_kwargs_from_configdict(CONFIG, doc=_GDPoptAlgorithm.solve.__doc__)
     def solve(self, model, **kwds):
         return super().solve(model, **kwds)
 
+    @staticmethod
+    def _quadratic_matrix(repn):
+        var_to_idx = ComponentMap()
+        for _coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars):
+            if v1 not in var_to_idx:
+                var_to_idx[v1] = len(var_to_idx)
+            if v2 not in var_to_idx:
+                var_to_idx[v2] = len(var_to_idx)
+
+        q_matrix = [
+            [0.0 for _ in range(len(var_to_idx))] for _ in range(len(var_to_idx))
+        ]
+        for coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars):
+            coef_val = value(coef, exception=False)
+            if coef_val is None:
+                return None
+            idx1 = var_to_idx[v1]
+            idx2 = var_to_idx[v2]
+            if v1 is v2:
+                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
+
+    @staticmethod
+    def _symmetric_matrix_eigenvalues(matrix):
+        if numpy_available:
+            return np.linalg.eigvalsh(matrix)
+
+        order = len(matrix)
+        if order <= 1:
+            return [matrix[0][0]] if order else []
+
+        matrix = [row[:] for row in matrix]
+        scale = max(abs(val) for row in matrix for val in row)
+        if scale == 0:
+            return [0.0 for _ in range(order)]
+
+        rotation_tolerance = 1e-12 * scale
+        max_rotations = max(1, 50 * order * order)
+        for _ in range(max_rotations):
+            pivot_i, pivot_j, offdiag = 0, 1, abs(matrix[0][1])
+            for i in range(order - 1):
+                for j in range(i + 1, order):
+                    candidate = abs(matrix[i][j])
+                    if candidate > offdiag:
+                        pivot_i, pivot_j, offdiag = i, j, candidate
+
+            if offdiag <= rotation_tolerance:
+                return [matrix[i][i] for i in range(order)]
+
+            pivot = matrix[pivot_i][pivot_j]
+            diag_i = matrix[pivot_i][pivot_i]
+            diag_j = matrix[pivot_j][pivot_j]
+            if diag_i == diag_j:
+                tangent = copysign(1.0, pivot)
+            else:
+                tau = (diag_j - diag_i) / (2.0 * pivot)
+                tangent = copysign(1.0, tau) / (abs(tau) + sqrt(1.0 + tau**2))
+            cosine = 1.0 / sqrt(1.0 + tangent**2)
+            sine = tangent * cosine
+
+            matrix[pivot_i][pivot_i] = diag_i - tangent * pivot
+            matrix[pivot_j][pivot_j] = diag_j + tangent * pivot
+            matrix[pivot_i][pivot_j] = matrix[pivot_j][pivot_i] = 0.0
+
+            for k in range(order):
+                if k in (pivot_i, pivot_j):
+                    continue
+                elem_i = matrix[k][pivot_i]
+                elem_j = matrix[k][pivot_j]
+                matrix[k][pivot_i] = matrix[pivot_i][k] = (
+                    cosine * elem_i - sine * elem_j
+                )
+                matrix[k][pivot_j] = matrix[pivot_j][k] = (
+                    sine * elem_i + cosine * elem_j
+                )
+
+        return None
+
+    @staticmethod
+    def _quadratic_curvature(expr):
+        repn = generate_standard_repn(expr, quadratic=True)
+        if not repn.quadratic_coefs:
+            return 0
+
+        q_matrix = GDP_LOA_Solver._quadratic_matrix(repn)
+        if q_matrix is None:
+            return None
+
+        eigenvalue_tolerance = 1e-10
+        eigenvalues = GDP_LOA_Solver._symmetric_matrix_eigenvalues(q_matrix)
+        if eigenvalues is None:
+            return None
+        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 _problem_may_have_nonrigorous_dual_bound(self, model):
+        for obj in model.component_data_objects(
+            Objective, active=True, descend_into=True
+        ):
+            degree = obj.expr.polynomial_degree()
+            if degree is None or degree not in (0, 1, 2):
+                return True
+            if degree == 2:
+                curvature = self._quadratic_curvature(obj.expr)
+                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)
+        ):
+            degree = constr.body.polynomial_degree()
+            if degree in (0, 1):
+                continue
+            if degree is None or degree not in (2,):
+                return True
+            if constr.equality:
+                return True
+            curvature = self._quadratic_curvature(constr.body)
+            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
+
     def _log_citation(self, config):
         config.logger.info("\n" + """- LOA algorithm:
         Türkay, M; Grossmann, IE.

pyomo/contrib/gdpopt/tests/test_gdpopt.py

@@ -12,6 +12,7 @@
 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
 
@@ -20,6 +21,9 @@
 from pyomo.common.collections import Bunch
 from pyomo.common.config import ConfigDict, ConfigValue
 from pyomo.common.fileutils import import_file, PYOMO_ROOT_DIR
+from pyomo.contrib.gdpopt.gloa import GDP_GLOA_Solver
+import pyomo.contrib.gdpopt.loa as loa_module
+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 +48,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 +83,108 @@
 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._nonrigorous_dual_bound_possible = 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_may_have_nonrigorous_dual_bound(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(GDP_LOA_Solver()._problem_may_have_nonrigorous_dual_bound(m))
+
+    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)
+
+        solver = GDP_LOA_Solver()
+        self.assertFalse(solver._problem_may_have_nonrigorous_dual_bound(convex))
+        self.assertTrue(solver._problem_may_have_nonrigorous_dual_bound(nonconvex))
+
+    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(GDP_LOA_Solver()._problem_may_have_nonrigorous_dual_bound(m))
+
+    def test_loa_certifies_psd_quadratic_cross_terms_without_numpy(self):
+        convex = ConcreteModel()
+        convex.x = Var(bounds=(-2, 2))
+        convex.y = Var(bounds=(-2, 2))
+        convex.c = Constraint(expr=convex.x**2 + convex.x * convex.y + convex.y**2 <= 4)
+        convex.obj = Objective(expr=convex.x)
+
+        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)
+
+        original_numpy_available = loa_module.numpy_available
+        loa_module.numpy_available = False
+        try:
+            self.assertFalse(
+                GDP_LOA_Solver()._problem_may_have_nonrigorous_dual_bound(convex)
+            )
+            self.assertTrue(
+                GDP_LOA_Solver()._problem_may_have_nonrigorous_dual_bound(nonconvex)
+            )
+        finally:
+            loa_module.numpy_available = original_numpy_available
+
+    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(
         SolverFactory(mip_solver).available(), "MIP solver not available"
     )