From e972e583a5dc2af79db5b8f43798a6d94c31c1fe Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Mon, 11 May 2026 18:48:39 -0400 Subject: [PATCH 1/3] Fix nonrigorous crossed bound results (#3947) --- pyomo/contrib/gdpopt/algorithm_base_class.py | 34 ++++++++- pyomo/contrib/gdpopt/loa.py | 52 +++++++++++++ pyomo/contrib/gdpopt/tests/test_gdpopt.py | 74 ++++++++++++++++++- pyomo/contrib/mindtpy/algorithm_base_class.py | 20 ++++- .../mindtpy/global_outer_approximation.py | 1 + .../mindtpy/tests/test_mindtpy_no_discrete.py | 36 ++++++++- 6 files changed, 212 insertions(+), 5 deletions(-) diff --git a/pyomo/contrib/gdpopt/algorithm_base_class.py b/pyomo/contrib/gdpopt/algorithm_base_class.py index 8463178e522..fac6e1b546a 100644 --- a/pyomo/contrib/gdpopt/algorithm_base_class.py +++ b/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 diff --git a/pyomo/contrib/gdpopt/loa.py b/pyomo/contrib/gdpopt/loa.py index 8e36bebfd22..b96e0623df1 100644 --- a/pyomo/contrib/gdpopt/loa.py +++ b/pyomo/contrib/gdpopt/loa.py @@ -74,12 +74,64 @@ 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_curvature(expr): + repn = generate_standard_repn(expr, quadratic=True) + curvature = 0 + for coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars): + if v1 is not v2: + return None + coef_val = value(coef, exception=False) + if coef_val is None: + return None + if coef_val > 0: + if curvature < 0: + return None + curvature = 1 + elif coef_val < 0: + if curvature > 0: + return None + curvature = -1 + return curvature + + 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. diff --git a/pyomo/contrib/gdpopt/tests/test_gdpopt.py b/pyomo/contrib/gdpopt/tests/test_gdpopt.py index 5a0fc30cf37..94dd39e5b81 100644 --- a/pyomo/contrib/gdpopt/tests/test_gdpopt.py +++ b/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,8 @@ 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 +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 +47,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 +82,74 @@ 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_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" ) diff --git a/pyomo/contrib/mindtpy/algorithm_base_class.py b/pyomo/contrib/mindtpy/algorithm_base_class.py index 11b34e47d8e..72c86ff4583 100644 --- a/pyomo/contrib/mindtpy/algorithm_base_class.py +++ b/pyomo/contrib/mindtpy/algorithm_base_class.py @@ -84,6 +84,8 @@ class _MindtPyAlgorithm: + _crossed_bounds_are_certified = True + def __init__(self, **kwds): """ This is a common init method for all the MindtPy algorithms, so that we @@ -104,6 +106,7 @@ def __init__(self, **kwds): self.timing = Bunch() self.curr_int_sol = [] self.should_terminate = False + self._nonrigorous_crossed_bounds = False self.integer_list = [] # Dictionary {integer solution (tuple): [cuts begin index, cuts end index] (list)} self.integer_solution_to_cuts_index = dict() @@ -2451,7 +2454,14 @@ def setup_regularization_main(self): ) def update_result(self): - if self.objective_sense == minimize: + if self._nonrigorous_crossed_bounds: + if self.objective_sense == minimize: + self.results.problem.lower_bound = float('-inf') + self.results.problem.upper_bound = self.primal_bound + else: + self.results.problem.lower_bound = self.primal_bound + self.results.problem.upper_bound = float('inf') + elif self.objective_sense == minimize: self.results.problem.lower_bound = self.dual_bound self.results.problem.upper_bound = self.primal_bound else: @@ -3032,6 +3042,7 @@ def solve(self, model, **kwds): kwds.pop('options', {}), preserve_implicit=True ) config.set_value(kwds) + self._nonrigorous_crossed_bounds = False self.set_up_logger() new_logging_level = logging.INFO if config.tee else None with lower_logger_level_to(config.logger, new_logging_level): @@ -3349,6 +3360,13 @@ def add_regularization(self): def bounds_converged(self): # Check bound convergence + if self.abs_gap < 0 and not self._crossed_bounds_are_certified: + self._nonrigorous_crossed_bounds = True + self.config.logger.info( + 'MindtPy exiting on crossed bounds without a certified dual bound.' + ) + self.results.solver.termination_condition = tc.feasible + return True if self.abs_gap <= self.config.absolute_bound_tolerance: self.config.logger.info( 'MindtPy exiting on bound convergence. ' diff --git a/pyomo/contrib/mindtpy/global_outer_approximation.py b/pyomo/contrib/mindtpy/global_outer_approximation.py index 9f3f70bde10..e8e39067d48 100644 --- a/pyomo/contrib/mindtpy/global_outer_approximation.py +++ b/pyomo/contrib/mindtpy/global_outer_approximation.py @@ -35,6 +35,7 @@ class MindtPy_GOA_Solver(_MindtPyAlgorithm): """ CONFIG = _get_MindtPy_GOA_config() + _crossed_bounds_are_certified = False def check_config(self): config = self.config diff --git a/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py b/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py index 1ba46b08d75..107250dfe37 100644 --- a/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py +++ b/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py @@ -8,9 +8,11 @@ # ____________________________________________________________________________________ +import logging from unittest.mock import MagicMock, patch -from pyomo.opt import TerminationCondition as tc, SolverStatus +from pyomo.common.collections import Bunch +from pyomo.opt import SolverResults, TerminationCondition as tc, SolverStatus import pyomo.common.unittest as unittest from pyomo.environ import ( @@ -264,6 +266,38 @@ def __init__(self, **kwargs): object.__setattr__(self, k, v) +class TestMindtPyGOAResults(unittest.TestCase): + def test_goa_crossed_bounds_are_not_reported_as_global_optimal(self): + from pyomo.contrib.mindtpy.global_outer_approximation import MindtPy_GOA_Solver + + solver = MindtPy_GOA_Solver() + solver.config = solver.CONFIG() + solver.config.absolute_bound_tolerance = 1e-6 + solver.config.relative_bound_tolerance = 1e-6 + solver.config.logger = logging.getLogger(__name__) + solver.results = SolverResults() + solver.objective_sense = maximize + solver.primal_bound = 1011.6577899409375 + solver.dual_bound = 1000.0 + solver.best_solution_found = object() + solver.update_gap() + + self.assertTrue(solver.bounds_converged()) + self.assertIs(solver.results.solver.termination_condition, tc.feasible) + + solver.timing = Bunch(total=0.0) + solver.mip_iter = 1 + solver.nlp_infeasible_counter = 0 + solver.best_solution_found_time = None + solver.primal_integral = 0.0 + solver.dual_integral = 0.0 + solver.primal_dual_gap_integral = 0.0 + solver.update_result() + + self.assertEqual(solver.results.problem.lower_bound, 1011.6577899409375) + self.assertEqual(solver.results.problem.upper_bound, float('inf')) + + class TestMirrorDirectSolveResults(unittest.TestCase): """Unit tests for _mirror_direct_solve_results covering all branches.""" From 3fa0b4bcceb012d34e22022d32ac10e27e2ac0fb Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Mon, 11 May 2026 19:43:00 -0400 Subject: [PATCH 2/3] Certify LOA convex quadratic cross terms (#3947) --- pyomo/contrib/gdpopt/loa.py | 38 +++++++++++++++++++++++ pyomo/contrib/gdpopt/tests/test_gdpopt.py | 9 ++++++ 2 files changed, 47 insertions(+) diff --git a/pyomo/contrib/gdpopt/loa.py b/pyomo/contrib/gdpopt/loa.py index b96e0623df1..03cf9c1161c 100644 --- a/pyomo/contrib/gdpopt/loa.py +++ b/pyomo/contrib/gdpopt/loa.py @@ -12,6 +12,7 @@ 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 ( @@ -84,6 +85,43 @@ def solve(self, model, **kwds): @staticmethod def _quadratic_curvature(expr): repn = generate_standard_repn(expr, quadratic=True) + if not repn.quadratic_coefs: + return 0 + + if numpy_available: + 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 = np.zeros((len(var_to_idx), 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 + + eigenvalue_tolerance = 1e-10 + eigenvalues = np.linalg.eigvalsh(q_matrix) + is_psd = not np.any(eigenvalues < -eigenvalue_tolerance) + is_nsd = not np.any(eigenvalues > eigenvalue_tolerance) + if is_psd and is_nsd: + return 0 + if is_psd: + return 1 + if is_nsd: + return -1 + return None + curvature = 0 for coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars): if v1 is not v2: diff --git a/pyomo/contrib/gdpopt/tests/test_gdpopt.py b/pyomo/contrib/gdpopt/tests/test_gdpopt.py index 94dd39e5b81..28594b4acb7 100644 --- a/pyomo/contrib/gdpopt/tests/test_gdpopt.py +++ b/pyomo/contrib/gdpopt/tests/test_gdpopt.py @@ -137,6 +137,15 @@ def test_loa_distinguishes_basic_quadratic_convexity(self): 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_gloa_crossed_bounds_preserve_certified_optimal_behavior(self): solver = GDP_GLOA_Solver() config = self._setup_crossed_bound_state(solver) From ce6b6d63c7ea7e9d06f08fdd190a485fab593980 Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Mon, 11 May 2026 21:57:20 -0400 Subject: [PATCH 3/3] Support LOA quadratic certification without NumPy (#3947) --- pyomo/contrib/gdpopt/loa.py | 148 +++++++++++++++------- pyomo/contrib/gdpopt/tests/test_gdpopt.py | 26 ++++ 2 files changed, 125 insertions(+), 49 deletions(-) diff --git a/pyomo/contrib/gdpopt/loa.py b/pyomo/contrib/gdpopt/loa.py index 03cf9c1161c..d15428550b3 100644 --- a/pyomo/contrib/gdpopt/loa.py +++ b/pyomo/contrib/gdpopt/loa.py @@ -8,7 +8,7 @@ # ____________________________________________________________________________________ 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 @@ -82,62 +82,112 @@ class GDP_LOA_Solver(_GDPoptAlgorithm, _OAAlgorithmMixIn): 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 - if numpy_available: - 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 = np.zeros((len(var_to_idx), 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 - - eigenvalue_tolerance = 1e-10 - eigenvalues = np.linalg.eigvalsh(q_matrix) - is_psd = not np.any(eigenvalues < -eigenvalue_tolerance) - is_nsd = not np.any(eigenvalues > eigenvalue_tolerance) - if is_psd and is_nsd: - return 0 - if is_psd: - return 1 - if is_nsd: - return -1 + q_matrix = GDP_LOA_Solver._quadratic_matrix(repn) + if q_matrix is None: return None - curvature = 0 - for coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars): - if v1 is not v2: - return None - coef_val = value(coef, exception=False) - if coef_val is None: - return None - if coef_val > 0: - if curvature < 0: - return None - curvature = 1 - elif coef_val < 0: - if curvature > 0: - return None - curvature = -1 - return curvature + 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( diff --git a/pyomo/contrib/gdpopt/tests/test_gdpopt.py b/pyomo/contrib/gdpopt/tests/test_gdpopt.py index 28594b4acb7..97f9959728f 100644 --- a/pyomo/contrib/gdpopt/tests/test_gdpopt.py +++ b/pyomo/contrib/gdpopt/tests/test_gdpopt.py @@ -22,6 +22,7 @@ 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, @@ -146,6 +147,31 @@ def test_loa_certifies_psd_quadratic_cross_terms(self): 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)