diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py index 643c2faf180..132b642d981 100644 --- a/pyomo/devel/initialization/examples/init_polynomial_ex.py +++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py @@ -55,6 +55,6 @@ def global_init_ex(): if __name__ == '__main__': # stat, x = lp_init_ex() - # stat, x = pwl_init_ex() - stat, x = global_init_ex() + stat, x = pwl_init_ex() + # stat, x = global_init_ex() print(stat, round(x, 4)) diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py index 25101b7d8d7..3743068c8f2 100644 --- a/pyomo/devel/initialization/global_init.py +++ b/pyomo/devel/initialization/global_init.py @@ -18,19 +18,27 @@ logger = logging.getLogger(__name__) -def _initialize_with_global_solver( - nlp: BlockData, global_solver: SolverBase, nlp_solver: SolverBase -): +def _set_global_solver_solution_limit(global_solver): + # Set solver specific option for solution limit if isinstance(global_solver, (ScipDirect, ScipPersistent)): - opts = {'limits/solutions': 1} + global_solver.config.solver_options['limits/solutions'] = 1 elif isinstance(global_solver, (GurobiDirectMINLP,)): - opts = {'SolutionLimit': 1} + global_solver.config.solver_options['SolutionLimit'] = 1 else: + # Raise error if solver is not currently implemented raise NotImplementedError( 'Currently, the initialization module only works with new solver ' 'interfaces, so the global solvers are limited to ScipDirect, ' 'ScipPersistent, and GurobiDirectMINLP.' ) + + +def _initialize_with_global_solver( + nlp: BlockData, global_solver: SolverBase, nlp_solver: SolverBase +): + # Set solution limit + _set_global_solver_solution_limit(global_solver) + # Check if time limit is provided for global solver if global_solver.config.time_limit is None: logger.warning( @@ -40,10 +48,7 @@ def _initialize_with_global_solver( ) res = global_solver.solve( - nlp, - load_solutions=False, - raise_exception_on_nonoptimal_result=False, - solver_options=opts, + nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False ) logger.info( f'solved NLP with {global_solver.name}: {res.solution_status}, {res.termination_condition}' diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py index 674e6547f38..0c12adadb96 100644 --- a/pyomo/devel/initialization/initialize.py +++ b/pyomo/devel/initialization/initialize.py @@ -97,6 +97,7 @@ def initialize_with_piecewise_linear_approximation( nlp_solver: SolverBase | None = None, mip_solver: SolverBase | None = None, default_bound: float = 1.0e8, + num_initial_points: int = 2, max_pwl_refinement_iter: int = 100, num_pwl_cons_to_refine_per_iter: int = 5, aggressive_substitution: bool = True, @@ -163,6 +164,7 @@ def initialize_with_piecewise_linear_approximation( mip_solver=mip_solver, nlp_solver=nlp_solver, default_bound=default_bound, + num_initial_points=num_initial_points, max_iter=max_pwl_refinement_iter, num_cons_to_refine_per_iter=num_pwl_cons_to_refine_per_iter, aggressive_substitution=aggressive_substitution, diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py index a8c98afb8db..c0239c11b03 100644 --- a/pyomo/devel/initialization/pwl_init.py +++ b/pyomo/devel/initialization/pwl_init.py @@ -19,6 +19,11 @@ PiecewiseLinearExpression, ) from pyomo.contrib.piecewise.piecewise_linear_function import PiecewiseLinearFunction +from pyomo.contrib.solver.solvers.scip.scip_direct import ScipDirect +from pyomo.contrib.solver.solvers.scip.scip_persistent import ScipPersistent +from pyomo.contrib.solver.solvers.gurobi.gurobi_direct_minlp import GurobiDirectMINLP +from pyomo.contrib.solver.solvers.gurobi.gurobi_persistent import GurobiPersistent +from pyomo.contrib.solver.solvers.highs import Highs from pyomo.contrib.solver.common.base import SolverBase from pyomo.contrib.solver.common.results import SolutionStatus from pyomo.core.base.block import BlockData @@ -230,11 +235,12 @@ def _refine_pwl_approx( violations.sort(key=lambda i: i[0], reverse=True) if len(violations) == 0: - raise RuntimeError( + logger.warning( 'We have not found a feasible solution to the problem yet, but the ' 'solution to piecewise linear approximation did not have any violations, ' - 'so there is nothing to refine.' + 'so there is nothing to refine. Ending refinement loop.' ) + return False tol = 1e-5 if math.isclose(violations[0][0], 0, abs_tol=tol): @@ -256,12 +262,31 @@ def _refine_pwl_approx( cons = pwl_expr_to_con_map.pop(e1) pwl_expr_to_con_map[e2] = cons + return True + + +def _set_mip_solver_solution_limit(mip_solver): + # Set solver specific option for solution limit + if isinstance(mip_solver, (ScipDirect, ScipPersistent)): + mip_solver.config.solver_options['limits/solutions'] = 1 + elif isinstance(mip_solver, (GurobiDirectMINLP, GurobiPersistent)): + mip_solver.config.solver_options['SolutionLimit'] = 1 + elif isinstance(mip_solver, Highs): + mip_solver.config.solver_options['mip_max_improving_sols'] = 1 + else: + raise NotImplementedError( + 'Currently, the initialization module only works with new solver ' + 'interfaces, so the mip solvers are limited to Highs, ScipDirect, ' + 'ScipPersistent, GurobiDirectMINLP, and GurobiPersistent.' + ) + def _initialize_with_piecewise_linear_approximation( nlp: BlockData, mip_solver: SolverBase, nlp_solver: SolverBase, default_bound=1.0e8, + num_initial_points=2, max_iter=100, num_cons_to_refine_per_iter=5, aggressive_substitution=True, @@ -297,7 +322,7 @@ def _initialize_with_piecewise_linear_approximation( # build the PWL approximation trans = pyo.TransformationFactory('contrib.piecewise.nonlinear_to_pwl') - trans.apply_to(pwl, num_points=2, additively_decompose=False) + trans.apply_to(pwl, num_points=num_initial_points, additively_decompose=False) logger.info('replaced nonlinear expressions with piecewise linear expressions') """ @@ -311,6 +336,7 @@ def _initialize_with_piecewise_linear_approximation( pwl_expr_to_con_map = _get_pwl_constraints(pwl) solved = False last_nlp_res = None + for _iter in range(max_iter): logger.info(f'PWL initialization: iter {_iter}') @@ -324,6 +350,8 @@ def _initialize_with_piecewise_linear_approximation( del _pwl.orig_vars logger.info('applied the disaggregated logarithmic transformation') + if max_iter == 1: + _set_mip_solver_solution_limit(mip_solver) # solve the MILP res = mip_solver.solve( _pwl, load_solutions=False, raise_exception_on_nonoptimal_result=False @@ -336,15 +364,6 @@ def _initialize_with_piecewise_linear_approximation( for ov, nv in zip(orig_vars, new_vars): ov.set_value(nv.value, skip_validation=True) - # refine the PWL approximation - _refine_pwl_approx( - pwl, - pwl_expr_to_con_map=pwl_expr_to_con_map, - num_to_refine=num_cons_to_refine_per_iter, - bounds_tol=bounds_tol, - ) - logger.info('refined PWL approximation') - # try solving the NLP res = nlp_solver.solve( nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False @@ -356,6 +375,23 @@ def _initialize_with_piecewise_linear_approximation( res.solution_loader.load_vars() break + # load the variable values back into orig_vars + for ov, nv in zip(orig_vars, new_vars): + ov.set_value(nv.value, skip_validation=True) + + # refine the PWL approximation, use refined check to decide whether to break + refined = _refine_pwl_approx( + pwl, + pwl_expr_to_con_map=pwl_expr_to_con_map, + num_to_refine=num_cons_to_refine_per_iter, + bounds_tol=bounds_tol, + ) + + if refined: + logger.info('refined PWL approximation') + else: + break + if not solved: logger.warning('initialization was not successful via PWL approximation') diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py index 042bb890828..d107aea809b 100644 --- a/pyomo/devel/initialization/tests/test_initialization.py +++ b/pyomo/devel/initialization/tests/test_initialization.py @@ -17,6 +17,8 @@ pwl_init_ex, global_init_ex, ) +from pyomo.devel.initialization.pwl_init import _set_mip_solver_solution_limit +from pyomo.devel.initialization.global_init import _set_global_solver_solution_limit from pyomo.common import unittest from pyomo.common.dependencies import scipy_available from pyomo.contrib.solver.common.factory import SolverFactory @@ -31,6 +33,7 @@ scip = SolverFactory('scip_direct') ipopt = SolverFactory('ipopt') highs = SolverFactory('highs') +gurobi = SolverFactory('gurobi_direct_minlp') class MockNLPSolver(SolverBase): @@ -209,6 +212,35 @@ def test_pwl_init(self): aggressive_substitution=False, ) + @unittest.skipUnless(highs.available(), 'highs is not available') + def test_pwl_init_single(self): + """ + Same as test_pwl_init but with single iteration + """ + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-15, 5)) + m.c = pyo.Constraint(expr=(m.x + 7) * (m.x + 5) * (m.x - 4) + 200 == 0) + m.obj = pyo.Objective(expr=m.x) + + # all the actual testing happens in the MockNLPSolver + nlp_solver = MockNLPSolver( + varlist=[m.x], + sol_map={ + 0: ([None], 0, 0), + 1: ([-9.920096055464825], 1e-4, 1e-4), + 2: ([-9.920096055464825], 1e-4, 1e-4), + }, + ) + mip_solver = SolverFactory('highs') + results = ini.initialize_with_piecewise_linear_approximation( + nlp=m, + nlp_solver=nlp_solver, + mip_solver=mip_solver, + num_initial_points=64, + max_pwl_refinement_iter=1, + aggressive_substitution=False, + ) + @unittest.skipUnless(highs.available(), 'highs is not available') @unittest.skipUnless(ipopt.available(), 'ipopt is not available') def test_pwl_ineq(self): @@ -227,6 +259,63 @@ def test_pwl_ineq(self): self.assertEqual(results.solution_status, SolutionStatus.optimal) self.assertAlmostEqual(results.incumbent_objective, 1, 5) + # Test solver support for global and pwl with no refinement + @unittest.skipUnless(highs.available(), 'highs is not available') + @unittest.skipUnless(scip.available(), 'scip is not available') + @unittest.skipUnless(gurobi.available(), 'gurobi is not available') + def test_global_solution_limit(self): + # Check supported solvers + solver_list = ["gurobi_direct_minlp", "scip_direct", "scip_persistent"] + + for solver in solver_list: + global_solver = SolverFactory(solver) + _set_global_solver_solution_limit(global_solver) + global_solver_opts = global_solver.config.solver_options + if solver == "gurobi_direct_minlp": + self.assertEqual(global_solver_opts["SolutionLimit"], 1) + elif solver in {"scip_direct", "scip_persistent"}: + self.assertEqual(global_solver_opts["limits/solutions"], 1) + # Check unsupported solver + with self.assertRaisesRegex( + NotImplementedError, + '.*Currently, the initialization module only works with new solver.*', + ): + wrong_solver = SolverFactory("highs") + _set_global_solver_solution_limit(wrong_solver) + + @unittest.skipUnless(highs.available(), 'highs is not available') + @unittest.skipUnless(scip.available(), 'scip is not available') + @unittest.skipUnless(gurobi.available(), 'gurobi is not available') + def test_pwl_solver_solution_limit(self): + + # Check supported solvers + solver_list = [ + "scip_direct", + "scip_persistent", + "gurobi_direct_minlp", + "gurobi_persistent", + "highs", + ] + + for solver in solver_list: + mip_solver = SolverFactory(solver) + _set_mip_solver_solution_limit(mip_solver) + mip_solver_opts = mip_solver.config.solver_options + if solver in {"gurobi_direct_minlp", "gurobi_persistent"}: + self.assertEqual(mip_solver_opts['SolutionLimit'], 1) + elif solver in {"scip_direct", "scip_persistent"}: + self.assertEqual(mip_solver_opts['limits/solutions'], 1) + elif solver == "highs": + self.assertEqual(mip_solver_opts['mip_max_improving_sols'], 1) + + # Check unsupported solver + with self.assertRaisesRegex( + NotImplementedError, + '.*Currently, the initialization module only works with new solver.*', + ): + wrong_solver = SolverFactory("multistart") + _set_mip_solver_solution_limit(wrong_solver) + if __name__ == '__main__': import logging