diff --git a/pyomo/contrib/gdpopt/algorithm_base_class.py b/pyomo/contrib/gdpopt/algorithm_base_class.py
index 8463178e522..a2bb8fcc7fa 100644
--- a/pyomo/contrib/gdpopt/algorithm_base_class.py
+++ b/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
diff --git a/pyomo/contrib/gdpopt/config_options.py b/pyomo/contrib/gdpopt/config_options.py
index 2e8e3ad309c..33eb01bd3d0 100644
--- a/pyomo/contrib/gdpopt/config_options.py
+++ b/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):
diff --git a/pyomo/contrib/gdpopt/convexity.py b/pyomo/contrib/gdpopt/convexity.py
new file mode 100644
index 00000000000..5eaef53331a
--- /dev/null
+++ b/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
diff --git a/pyomo/contrib/gdpopt/loa.py b/pyomo/contrib/gdpopt/loa.py
index 3f9db839b93..151caf91ccc 100644
--- a/pyomo/contrib/gdpopt/loa.py
+++ b/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__)
diff --git a/pyomo/contrib/gdpopt/tests/test_gdpopt.py b/pyomo/contrib/gdpopt/tests/test_gdpopt.py
index 95dc57f763a..b22c3d5b77e 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
 
@@ -19,7 +20,11 @@
 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"
     )
diff --git a/pyomo/contrib/mindtpy/algorithm_base_class.py b/pyomo/contrib/mindtpy/algorithm_base_class.py
index 11b34e47d8e..128e51d04c3 100644
--- a/pyomo/contrib/mindtpy/algorithm_base_class.py
+++ b/pyomo/contrib/mindtpy/algorithm_base_class.py
@@ -55,6 +55,7 @@
 from pyomo.contrib.gdpopt.solve_discrete_problem import (
     distinguish_mip_infeasible_or_unbounded,
 )
+from pyomo.contrib.gdpopt.convexity import model_is_not_certified_convex
 from pyomo.contrib.mindtpy.util import (
     generate_norm1_objective_function,
     generate_norm2sq_objective_function,
@@ -84,6 +85,10 @@
 
 
 class _MindtPyAlgorithm:
+    # OA and ECP need a convex model for their linearizations to provide a valid
+    # relaxation. GOA uses McCormick envelopes and does not have this requirement.
+    _requires_model_convexity = False
+
     def __init__(self, **kwds):
         """
         This is a common init method for all the MindtPy algorithms, so that we
@@ -104,6 +109,8 @@ def __init__(self, **kwds):
         self.timing = Bunch()
         self.curr_int_sol = []
         self.should_terminate = False
+        self._bounds_crossed_without_certified_convexity = False
+        self._model_is_not_certified_convex = False
         self.integer_list = []
         # Dictionary {integer solution (tuple): [cuts begin index, cuts end index] (list)}
         self.integer_solution_to_cuts_index = dict()
@@ -2451,7 +2458,14 @@ def setup_regularization_main(self):
             )
 
     def update_result(self):
-        if self.objective_sense == minimize:
+        if self._bounds_crossed_without_certified_convexity:
+            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 +3046,8 @@ def solve(self, model, **kwds):
             kwds.pop('options', {}), preserve_implicit=True
         )
         config.set_value(kwds)
+        self._bounds_crossed_without_certified_convexity = False
+        self._model_is_not_certified_convex = 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):
@@ -3059,6 +3075,12 @@ def solve(self, model, **kwds):
                 self.results.problem.number_of_objectives = (
                     self._original_model_num_active_objectives
                 )
+                self._model_is_not_certified_convex = (
+                    self._requires_model_convexity
+                    and model_is_not_certified_convex(
+                        self.original_model, config.eigenvalue_tolerance
+                    )
+                )
 
                 # Validate the model to ensure that MindtPy is able to solve it.
                 if not self.model_is_valid():
@@ -3349,6 +3371,17 @@ def add_regularization(self):
 
     def bounds_converged(self):
         # Check bound convergence
+        if (
+            self.abs_gap < 0
+            and self._requires_model_convexity
+            and self._model_is_not_certified_convex
+        ):
+            self._bounds_crossed_without_certified_convexity = True
+            self.config.logger.info(
+                'MindtPy exiting on crossed bounds without certified model convexity.'
+            )
+            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/config_options.py b/pyomo/contrib/mindtpy/config_options.py
index 0ac163ccfcc..8391e3fc25b 100644
--- a/pyomo/contrib/mindtpy/config_options.py
+++ b/pyomo/contrib/mindtpy/config_options.py
@@ -13,6 +13,7 @@
     ConfigBlock,
     ConfigValue,
     In,
+    NonNegativeFloat,
     PositiveFloat,
     PositiveInt,
     NonNegativeInt,
@@ -674,6 +675,14 @@ def _add_tolerance_configs(CONFIG):
             ':math:`|Primal Bound - Dual Bound| / (1e-10 + |Primal Bound|) <= relative tolerance`',
         ),
     )
+    CONFIG.declare(
+        'eigenvalue_tolerance',
+        ConfigValue(
+            default=1e-10,
+            domain=NonNegativeFloat,
+            description='Numerical tolerance for eigenvalue-based convexity checks.',
+        ),
+    )
     CONFIG.declare(
         'small_dual_tolerance',
         ConfigValue(
diff --git a/pyomo/contrib/mindtpy/extended_cutting_plane.py b/pyomo/contrib/mindtpy/extended_cutting_plane.py
index a209125fcd0..2b7de4fbec3 100644
--- a/pyomo/contrib/mindtpy/extended_cutting_plane.py
+++ b/pyomo/contrib/mindtpy/extended_cutting_plane.py
@@ -35,6 +35,8 @@ class MindtPy_ECP_Solver(_MindtPyAlgorithm):
     """
 
     CONFIG = _get_MindtPy_ECP_config()
+    # ECP linearizations only provide a valid relaxation for convex models.
+    _requires_model_convexity = True
 
     def MindtPy_iteration_loop(self):
         """Main loop for MindtPy Algorithms.
diff --git a/pyomo/contrib/mindtpy/global_outer_approximation.py b/pyomo/contrib/mindtpy/global_outer_approximation.py
index 9f3f70bde10..93569acdb55 100644
--- a/pyomo/contrib/mindtpy/global_outer_approximation.py
+++ b/pyomo/contrib/mindtpy/global_outer_approximation.py
@@ -35,6 +35,8 @@ class MindtPy_GOA_Solver(_MindtPyAlgorithm):
     """
 
     CONFIG = _get_MindtPy_GOA_config()
+    # McCormick envelopes provide a valid relaxation without requiring the
+    # original model to be convex.
 
     def check_config(self):
         config = self.config
diff --git a/pyomo/contrib/mindtpy/outer_approximation.py b/pyomo/contrib/mindtpy/outer_approximation.py
index cdd54b4a1f5..f3575f82c63 100644
--- a/pyomo/contrib/mindtpy/outer_approximation.py
+++ b/pyomo/contrib/mindtpy/outer_approximation.py
@@ -36,6 +36,8 @@ class MindtPy_OA_Solver(_MindtPyAlgorithm):
     """
 
     CONFIG = _get_MindtPy_OA_config()
+    # OA linearizations only provide a valid relaxation for convex models.
+    _requires_model_convexity = True
 
     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..72aa461d5da 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 (
@@ -25,6 +27,10 @@
     maximize,
     value,
 )
+from pyomo.contrib.mindtpy.algorithm_base_class import _MindtPyAlgorithm
+from pyomo.contrib.mindtpy.extended_cutting_plane import MindtPy_ECP_Solver
+from pyomo.contrib.mindtpy.global_outer_approximation import MindtPy_GOA_Solver
+from pyomo.contrib.mindtpy.outer_approximation import MindtPy_OA_Solver
 
 required_nlp_solvers = 'ipopt'
 # Open-source (or generally available) solver pair used by MindtPy tests that
@@ -264,6 +270,78 @@ def __init__(self, **kwargs):
             object.__setattr__(self, k, v)
 
 
+def _crossed_bound_solver(solver, model_is_not_certified_convex):
+    """Set up a solver instance with maximization bounds that have crossed."""
+    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._model_is_not_certified_convex = model_is_not_certified_convex
+    solver.update_gap()
+    return solver
+
+
+def _finalize_result(solver):
+    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()
+
+
+class TestMindtPyCrossedBoundResults(unittest.TestCase):
+    def test_oa_crossed_bounds_are_not_reported_as_global_optimal(self):
+        solver = _crossed_bound_solver(MindtPy_OA_Solver(), True)
+
+        self.assertTrue(solver.bounds_converged())
+        self.assertIs(solver.results.solver.termination_condition, tc.feasible)
+
+        _finalize_result(solver)
+
+        self.assertEqual(solver.results.problem.lower_bound, 1011.6577899409375)
+        self.assertEqual(solver.results.problem.upper_bound, float('inf'))
+
+    def test_ecp_crossed_bounds_are_not_reported_as_global_optimal(self):
+        solver = _crossed_bound_solver(MindtPy_ECP_Solver(), True)
+
+        self.assertTrue(solver.bounds_converged())
+        self.assertIs(solver.results.solver.termination_condition, tc.feasible)
+
+    def test_oa_crossed_bounds_on_convex_model_stay_optimal(self):
+        """A convex model gives OA a rigorous dual bound, so crossing is tolerance."""
+        solver = _crossed_bound_solver(MindtPy_OA_Solver(), False)
+
+        self.assertTrue(solver.bounds_converged())
+        self.assertIs(solver.results.solver.termination_condition, tc.optimal)
+
+        _finalize_result(solver)
+
+        # The certified path reports both bounds unchanged, rather than replacing
+        # the uncertified side with an infinite bound.
+        self.assertEqual(solver.results.problem.lower_bound, 1011.6577899409375)
+        self.assertEqual(solver.results.problem.upper_bound, 1000.0)
+
+    def test_goa_crossed_bounds_preserve_certified_optimal_behavior(self):
+        """GOA relaxes with McCormick envelopes, so its dual bound is rigorous."""
+        solver = _crossed_bound_solver(MindtPy_GOA_Solver(), True)
+
+        self.assertTrue(solver.bounds_converged())
+        self.assertIs(solver.results.solver.termination_condition, tc.optimal)
+
+    def test_algorithm_convexity_requirements(self):
+        self.assertTrue(MindtPy_OA_Solver._requires_model_convexity)
+        self.assertTrue(MindtPy_ECP_Solver._requires_model_convexity)
+        self.assertFalse(MindtPy_GOA_Solver._requires_model_convexity)
+
+
 class TestMirrorDirectSolveResults(unittest.TestCase):
     """Unit tests for _mirror_direct_solve_results covering all branches."""
 
@@ -272,8 +350,6 @@ class TestMirrorDirectSolveResults(unittest.TestCase):
     def _make_algorithm_stub(self):
         """Create a minimal stub of _MindtPyAlgorithm with only the fields
         needed by _mirror_direct_solve_results."""
-        from pyomo.contrib.mindtpy.algorithm_base_class import _MindtPyAlgorithm
-
         stub = MagicMock(spec=_MindtPyAlgorithm)
         stub.results = MagicMock()
         stub.results.solver = MagicMock()
@@ -506,8 +582,6 @@ def _make_algorithm(
         mip_constraint_polynomial_degree=None,
         mip_objective_polynomial_degree=None,
     ):
-        from pyomo.contrib.mindtpy.algorithm_base_class import _MindtPyAlgorithm
-
         algo = _MindtPyAlgorithm()
         algo.config = _SimpleNamespace(
             logger=MagicMock(),