diff --git a/pyomo/contrib/mindtpy/algorithm_base_class.py b/pyomo/contrib/mindtpy/algorithm_base_class.py index cb1c34e3054..6aac9a0f19d 100644 --- a/pyomo/contrib/mindtpy/algorithm_base_class.py +++ b/pyomo/contrib/mindtpy/algorithm_base_class.py @@ -252,13 +252,14 @@ def model_is_valid(self): This method performs a structural check on the working model. It determines if the problem is a true Mixed-Integer program. If no discrete variables are present, it serves as a short-circuit. - In short-circuit cases, the problem is solved immediately as an LP or NLP. + In short-circuit cases, the problem is solved immediately with the + configured MIP or NLP subsolver. Returns ------- bool True if the model has discrete variables and requires MindtPy iteration. - False if the model is purely continuous (LP or NLP). + False if the model is purely continuous (LP, QP, QCP, or NLP). Notes ----- @@ -270,19 +271,20 @@ def model_is_valid(self): This indicates the model is a valid MINLP for decomposition. 2. Continuous Model Handling (The "False" cases) - If the discrete variable list is empty, the model is "invalid" for MINLP. - The method then differentiates between LP and NLP structures. + If the discrete variable list is empty, the model is "invalid" for + MINLP. The method then classifies the continuous model as LP, QP, QCP, + or NLP and routes it directly to the configured MIP or NLP subsolver. 3. NLP Branch - The code checks the ``polynomial_degree`` of constraints and objectives. - If any degree is non-linear (not in ``mip_constraint_polynomial_degree``), - it is treated as a standard Nonlinear Program (NLP). - The ``config.nlp_solver`` is called to solve the original model directly. + If the model is structurally nonlinear beyond QP/QCP, or if the chosen + MIP solver cannot handle the required quadratic structure, the + ``config.nlp_solver`` is called to solve the original model directly. - 4. LP Branch - If all components are linear, it is treated as a Linear Program (LP). - The ``config.mip_solver`` is utilized for the solution process. - Solutions are loaded directly back into the ``original_model``. + 4. MIP Branch + If the continuous model is LP, QP, or QCP and the chosen MIP solver + supports the required structure, ``config.mip_solver`` is used for the + direct solve. Solutions are loaded directly back into the + ``original_model``. In both continuous cases, the method returns False to bypass the main loop. This ensures MindtPy does not attempt decomposition on trivial continuous models. @@ -291,25 +293,53 @@ def model_is_valid(self): MindtPy = m.MindtPy_utils config = self.config - # Handle LP/NLP being passed to the solver + # Handle purely continuous models by short-circuiting to a direct solve prob = self.results.problem if len(MindtPy.discrete_variable_list) == 0: config.logger.info('Problem has no discrete decisions.') - obj = next(m.component_data_objects(ctype=Objective, active=True)) - if ( - any( - c.body.polynomial_degree() - not in self.mip_constraint_polynomial_degree - for c in MindtPy.constraint_list - ) - or obj.expr.polynomial_degree() - not in self.mip_objective_polynomial_degree - ): - config.logger.info( - 'Your model is a NLP (nonlinear program). ' - 'Using NLP solver %s to solve.' % config.nlp_solver + original_obj = next( + self.original_model.component_data_objects(ctype=Objective, active=True) + ) + problem_type_names = { + 'LP': 'LP (linear program)', + 'QP': 'QP (quadratic program)', + 'QCP': 'QCP (quadratically constrained program)', + 'NLP': 'NLP (nonlinear program)', + } + problem_type_articles = {'LP': 'an', 'QP': 'a', 'QCP': 'a', 'NLP': 'an'} + obj_degree = MindtPy.objective_polynomial_degree + if obj_degree is None: + working_obj = next( + m.component_data_objects(ctype=Objective, active=True) ) + obj_degree = working_obj.expr.polynomial_degree() + problem_type, solver_to_use, unsupported_structure = ( + self._classify_short_circuit_problem(MindtPy, obj_degree) + ) + if solver_to_use == 'nlp': + if problem_type == 'NLP': + config.logger.info( + 'Your model is %s %s. ' + 'Using NLP solver %s to solve.' + % ( + problem_type_articles[problem_type], + problem_type_names[problem_type], + config.nlp_solver, + ) + ) + else: + config.logger.info( + 'Your model is %s %s, but MIP solver %s does not support %s. ' + 'Using NLP solver %s to solve.' + % ( + problem_type_articles[problem_type], + problem_type_names[problem_type], + config.mip_solver, + unsupported_structure, + config.nlp_solver, + ) + ) update_solver_timelimit( self.nlp_opt, config.nlp_solver, self.timing, config ) @@ -322,12 +352,18 @@ def model_is_valid(self): if len(results.solution) > 0: self.original_model.solutions.load_from(results) - self._mirror_direct_solve_results(results=results, obj=obj, prob=prob) + self._mirror_direct_solve_results( + results=results, obj=original_obj, prob=prob + ) return False else: config.logger.info( - 'Your model is an LP (linear program). ' - 'Using LP solver %s to solve.' % config.mip_solver + 'Your model is %s %s. Using MIP solver %s to solve.' + % ( + problem_type_articles[problem_type], + problem_type_names[problem_type], + config.mip_solver, + ) ) if isinstance(self.mip_opt, PersistentSolver): self.mip_opt.set_instance(self.original_model) @@ -343,7 +379,9 @@ def model_is_valid(self): if len(results.solution) > 0: self.original_model.solutions.load_from(results) - self._mirror_direct_solve_results(results=results, obj=obj, prob=prob) + self._mirror_direct_solve_results( + results=results, obj=original_obj, prob=prob + ) return False # Set up dual value reporting @@ -359,12 +397,88 @@ def model_is_valid(self): # need to do some kind of transformation (Glover?) or throw an error message return True + def _classify_short_circuit_problem(self, MindtPy, obj_degree): + """Classify a no-discrete model for short-circuit direct solves. + + This classification is intentionally independent of + ``quadratic_strategy``. In the no-discrete short-circuit path we are + selecting a direct subsolver for the original model, not building the + MindtPy main problem. + + Parameters + ---------- + MindtPy : Block + The utility block attached to the working model + (``model.MindtPy_utils``). Must have ``has_quadratic_constraints`` + and ``has_nonquadratic_constraints`` boolean attributes set by + ``build_ordered_component_lists``. + obj_degree : int or None + Polynomial degree of the active objective expression, or ``None`` + for nonlinear (non-polynomial) objectives. + + Returns + ------- + tuple + ``(problem_type, solver_to_use, unsupported_structure)`` where + ``problem_type`` is one of ``LP``, ``QP``, ``QCP``, or ``NLP``, + ``solver_to_use`` is ``mip`` or ``nlp``, and + ``unsupported_structure`` is ``None`` unless a quadratic problem + must be routed to the NLP solver because the chosen MIP solver does + not support the required structure. + """ + has_quadratic_constraints = MindtPy.has_quadratic_constraints + has_nonquadratic_constraints = MindtPy.has_nonquadratic_constraints + + if has_nonquadratic_constraints or obj_degree not in (0, 1, 2): + return 'NLP', 'nlp', None + + if has_quadratic_constraints: + if not self._mip_solver_supports_capability('quadratic_constraint'): + return 'QCP', 'nlp', 'quadratic constraints' + if obj_degree == 2 and not self._mip_solver_supports_capability( + 'quadratic_objective' + ): + return 'QCP', 'nlp', 'quadratic objectives' + return 'QCP', 'mip', None + + if obj_degree == 2: + if not self._mip_solver_supports_capability('quadratic_objective'): + return 'QP', 'nlp', 'quadratic objectives' + return 'QP', 'mip', None + + return 'LP', 'mip', None + + def _mip_solver_supports_capability(self, capability): + """Return whether the selected MIP solver supports a model feature.""" + appsi_capabilities = { + 'appsi_cplex': {'quadratic_objective': True, 'quadratic_constraint': True}, + 'appsi_gurobi': {'quadratic_objective': True, 'quadratic_constraint': True}, + 'appsi_highs': { + # TODO: revisit if/when the APPSI HiGHS wrapper adds QP support. + 'quadratic_objective': False, + 'quadratic_constraint': False, + }, + } + solver_name = self.config.mip_solver + if solver_name in appsi_capabilities: + return appsi_capabilities[solver_name][capability] + + # APPSI solvers do not expose has_capability (that is a legacy-interface + # method). Any APPSI solver not in the dict above therefore reaches this + # branch and is treated conservatively: all quadratic features are assumed + # unsupported, so QP/QCP models are routed to the NLP solver. To enable + # MIP-path routing for a new APPSI solver, add it to appsi_capabilities. + has_capability = getattr(self.mip_opt, 'has_capability', None) + if has_capability is None: + return False + return has_capability(capability) + def _mirror_direct_solve_results(self, results, obj, prob): - """Mirror a direct (LP/NLP) solve result into MindtPy's results object. + """Mirror a direct short-circuit solve result into MindtPy results. This is used by `model_is_valid()` when the instance is purely continuous - (no discrete variables) and MindtPy short-circuits to a direct LP/NLP - solve. + (no discrete variables) and MindtPy short-circuits to a direct MIP or + NLP solve for an LP, QP, QCP, or NLP model. Parameters ---------- @@ -431,21 +545,31 @@ def build_ordered_component_lists(self, model): ) else: util_block.grey_box_list = [] - util_block.linear_constraint_list = list( - c - for c in util_block.constraint_list - if c.body.polynomial_degree() in self.mip_constraint_polynomial_degree - ) - util_block.nonlinear_constraint_list = list( - c - for c in util_block.constraint_list - if c.body.polynomial_degree() not in self.mip_constraint_polynomial_degree - ) + util_block.linear_constraint_list = [] + util_block.nonlinear_constraint_list = [] + util_block.has_quadratic_constraints = False + util_block.has_nonquadratic_constraints = False + for c in util_block.constraint_list: + degree = c.body.polynomial_degree() + if degree in self.mip_constraint_polynomial_degree: + util_block.linear_constraint_list.append(c) + else: + util_block.nonlinear_constraint_list.append(c) + + if degree == 2: + util_block.has_quadratic_constraints = True + elif degree not in (0, 1): + util_block.has_nonquadratic_constraints = True util_block.objective_list = list( model.component_data_objects( ctype=Objective, active=True, descend_into=(Block) ) ) + util_block.objective_polynomial_degree = None + if len(util_block.objective_list) == 1: + util_block.objective_polynomial_degree = util_block.objective_list[ + 0 + ].expr.polynomial_degree() # Identify the non-fixed variables in (potentially) active constraints and # objective functions diff --git a/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py b/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py index a9219ceb8b5..fed46e33122 100644 --- a/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py +++ b/pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py @@ -8,7 +8,7 @@ # ____________________________________________________________________________________ -from unittest.mock import MagicMock +from unittest.mock import MagicMock, patch from pyomo.opt import TerminationCondition as tc, SolverStatus import pyomo.common.unittest as unittest @@ -405,3 +405,316 @@ def test_solver_status_and_message_mirrored(self): self.assertEqual(algo.results.solver.status, SolverStatus.ok) self.assertEqual(algo.results.solver.termination_condition, tc.optimal) self.assertEqual(algo.results.solver.message, "All good") + + +class _FakeLegacyMIPSolver: + def __init__( + self, + *, + quadratic_objective=False, + quadratic_constraint=False, + termination_condition=tc.optimal, + ): + self.options = {} + self._capabilities = { + 'quadratic_objective': quadratic_objective, + 'quadratic_constraint': quadratic_constraint, + } + self.solve = MagicMock( + return_value=_make_direct_solve_results(termination_condition) + ) + + def has_capability(self, capability): + return self._capabilities[capability] + + +class _FakeSolver: + def __init__(self, termination_condition=tc.optimal): + self.options = {} + self.solve = MagicMock( + return_value=_make_direct_solve_results(termination_condition) + ) + + +def _make_direct_solve_results(termination_condition): + return _SimpleNamespace( + solution=[], + solver=_SimpleNamespace( + termination_condition=termination_condition, + status=SolverStatus.ok, + message=None, + ), + problem=_SimpleNamespace(lower_bound=None, upper_bound=None), + ) + + +class TestMindtPyShortCircuitRouting(unittest.TestCase): + def _make_algorithm( + self, + model, + mip_solver_name, + mip_opt, + nlp_opt=None, + 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(), + mip_solver=mip_solver_name, + nlp_solver='ipopt', + mip_solver_tee=False, + nlp_solver_tee=False, + mip_solver_args={}, + nlp_solver_args={}, + feasibility_norm='L1', + add_slack=False, + max_slack=0, + ) + algo.mip_constraint_polynomial_degree = ( + {0, 1} + if mip_constraint_polynomial_degree is None + else mip_constraint_polynomial_degree + ) + algo.mip_objective_polynomial_degree = ( + {0, 1} + if mip_objective_polynomial_degree is None + else mip_objective_polynomial_degree + ) + algo.original_model = model + algo.working_model = model.clone() + algo.create_utility_block(algo.working_model, 'MindtPy_utils') + algo.mip_opt = mip_opt + algo.nlp_opt = nlp_opt if nlp_opt is not None else _FakeSolver() + algo.mip_load_solutions = False + algo.nlp_load_solutions = False + algo._mirror_direct_solve_results = MagicMock() + algo.timing = _SimpleNamespace(main_timer_start_time=0.0) + return algo + + def _make_lp_model(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.c = Constraint(expr=m.x >= 1) + m.obj = Objective(expr=m.x, sense=minimize) + return m + + def _make_qp_model(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.c = Constraint(expr=m.x >= 1) + m.obj = Objective(expr=(m.x - 2) ** 2, sense=minimize) + return m + + def _make_qcp_model(self): + m = ConcreteModel() + m.x = Var(bounds=(0, 10)) + m.c = Constraint(expr=m.x**2 <= 4) + m.obj = Objective(expr=m.x, sense=maximize) + return m + + def _make_qcp_with_quadratic_objective_model(self): + m = ConcreteModel() + m.x = Var(bounds=(0, 10)) + m.c = Constraint(expr=m.x**2 <= 9) + m.obj = Objective(expr=(m.x - 1) ** 2, sense=minimize) + return m + + def _make_nlp_model(self): + m = ConcreteModel() + m.x = Var(bounds=(0, 10)) + m.c = Constraint(expr=m.x >= 1) + m.obj = Objective(expr=m.x**3, sense=minimize) + return m + + def test_short_circuit_lp_routes_to_mip(self): + algo = self._make_algorithm( + self._make_lp_model(), + mip_solver_name='glpk', + mip_opt=_FakeLegacyMIPSolver(), + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_called_once() + algo.nlp_opt.solve.assert_not_called() + + def test_short_circuit_qp_uses_legacy_mip_solver_when_supported(self): + algo = self._make_algorithm( + self._make_qp_model(), + mip_solver_name='gurobi', + mip_opt=_FakeLegacyMIPSolver(quadratic_objective=True), + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_called_once() + algo.nlp_opt.solve.assert_not_called() + self.assertIs( + algo._mirror_direct_solve_results.call_args.kwargs['obj'], + algo.original_model.obj, + ) + + def test_short_circuit_qp_uses_nlp_when_appsi_highs_lacks_quadratic_support(self): + algo = self._make_algorithm( + self._make_qp_model(), mip_solver_name='appsi_highs', mip_opt=_FakeSolver() + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_not_called() + algo.nlp_opt.solve.assert_called_once() + + def test_short_circuit_qcp_uses_appsi_cplex_when_supported(self): + algo = self._make_algorithm( + self._make_qcp_model(), mip_solver_name='appsi_cplex', mip_opt=_FakeSolver() + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_called_once() + algo.nlp_opt.solve.assert_not_called() + + def test_short_circuit_qcp_uses_nlp_when_legacy_mip_solver_lacks_qcp_support(self): + algo = self._make_algorithm( + self._make_qcp_model(), + mip_solver_name='cbc', + mip_opt=_FakeLegacyMIPSolver( + quadratic_objective=False, quadratic_constraint=False + ), + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_not_called() + algo.nlp_opt.solve.assert_called_once() + + def test_short_circuit_qcp_uses_nlp_when_solver_lacks_quadratic_objective(self): + algo = self._make_algorithm( + self._make_qcp_with_quadratic_objective_model(), + mip_solver_name='cbc', + mip_opt=_FakeLegacyMIPSolver( + quadratic_objective=False, quadratic_constraint=True + ), + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_not_called() + algo.nlp_opt.solve.assert_called_once() + + def test_short_circuit_nlp_uses_nlp_even_with_quadratic_capable_mip_solver(self): + algo = self._make_algorithm( + self._make_nlp_model(), + mip_solver_name='gurobi', + mip_opt=_FakeLegacyMIPSolver( + quadratic_objective=True, quadratic_constraint=True + ), + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_not_called() + algo.nlp_opt.solve.assert_called_once() + + def test_short_circuit_mip_failure_does_not_fallback_to_nlp(self): + algo = self._make_algorithm( + self._make_qp_model(), + mip_solver_name='gurobi', + mip_opt=_FakeLegacyMIPSolver( + quadratic_objective=True, termination_condition=tc.error + ), + ) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_called_once() + algo.nlp_opt.solve.assert_not_called() + + def test_short_circuit_qcp_detection_ignores_quadratic_strategy(self): + algo = self._make_algorithm( + self._make_qcp_model(), + mip_solver_name='cbc', + mip_opt=_FakeLegacyMIPSolver( + quadratic_objective=False, quadratic_constraint=False + ), + mip_constraint_polynomial_degree={0, 1, 2}, + mip_objective_polynomial_degree={0, 1, 2}, + ) + + util = algo.working_model.MindtPy_utils + self.assertEqual(len(util.linear_constraint_list), 1) + self.assertEqual(len(util.nonlinear_constraint_list), 0) + self.assertTrue(util.has_quadratic_constraints) + self.assertFalse(util.has_nonquadratic_constraints) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_not_called() + algo.nlp_opt.solve.assert_called_once() + + def _make_mixed_degree_model(self): + """Model with one quadratic and one cubic constraint (mixed degree).""" + m = ConcreteModel() + m.x = Var(bounds=(0, 10)) + m.c_quad = Constraint(expr=m.x**2 <= 9) + m.c_cubic = Constraint(expr=m.x**3 <= 27) + m.obj = Objective(expr=m.x, sense=minimize) + return m + + def test_short_circuit_mixed_degree_model_routes_to_nlp(self): + """A model with both quadratic and cubic constraints must route to NLP. + + has_quadratic_constraints and has_nonquadratic_constraints are both + True. The NLP short-circuit path should be taken regardless of whether + the configured MIP solver claims quadratic support, because the cubic + constraint makes this an NLP. + """ + algo = self._make_algorithm( + self._make_mixed_degree_model(), + mip_solver_name='gurobi', + mip_opt=_FakeLegacyMIPSolver( + quadratic_objective=True, quadratic_constraint=True + ), + ) + + util = algo.working_model.MindtPy_utils + self.assertTrue(util.has_quadratic_constraints) + self.assertTrue(util.has_nonquadratic_constraints) + + with patch( + 'pyomo.contrib.mindtpy.algorithm_base_class.update_solver_timelimit' + ): + self.assertFalse(algo.model_is_valid()) + + algo.mip_opt.solve.assert_not_called() + algo.nlp_opt.solve.assert_called_once()