Extend pwl approximation init method to support higher initial points and single solve.

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))

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}'

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,

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')
 

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