Update multistart solver, add method for multistart NLP initialization

pyomo/contrib/multistart/reinit.py

@@ -11,35 +11,49 @@
 
 import logging
 import random
-
-from pyomo.core import Var
+from pyomo.common.dependencies import numpy as np
+from pyomo.common.dependencies.scipy import stats
+from pyomo.util.vars_from_expressions import get_vars_from_components
 
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
-def rand(val, lb, ub):
-    return random.uniform(lb, ub)  # uniform distribution between lb and ub
+def rand(val, lb, ub, sampler):
+    # sample = sampler.rng.uniform(lb, ub)
+    sample = sampler.sample_scalar(lb, ub)  # uniform distribution between lb and ub
+    return sample
+
+
+def rand_vector(lbs, ubs, sampler):
+    lowers = lbs
+    uppers = ubs
+    # Generate vector of samples using sampler
+    samples = sampler.sample_vector(lowers, uppers)
+    return samples
 
 
-def midpoint_guess_and_bound(val, lb, ub):
+def midpoint_guess_and_bound(val, lb, ub, sampler=None):
     """Midpoint between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
     return (far_bound + val) / 2
 
 
-def rand_guess_and_bound(val, lb, ub):
+def rand_guess_and_bound(val, lb, ub, sampler):
     """Random choice between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
-    return random.uniform(val, far_bound)
+    if far_bound == ub:
+        return sampler.sample_scalar(val, far_bound)
+    else:
+        return sampler.sample_scalar(far_bound, val)
 
 
-def rand_distributed(val, lb, ub, divisions=9):
+def rand_distributed(val, lb, ub, sampler, divisions=9):
     """Random choice among evenly distributed set of values between bounds."""
     set_distributed_vals = linspace(lb, ub, divisions)
-    return random.choice(set_distributed_vals)
+    return sampler.rng.choice(set_distributed_vals)
 
 
-def simple_midpoint(val, lb, ub):
+def simple_midpoint(val, lb, ub, sampler=None):
     return (lb + ub) * 0.5
 
 
@@ -50,20 +64,23 @@ def linspace(lower, upper, n):
 
 strategies = {
     "rand": rand,
+    "rand_vector": rand_vector,
     "midpoint_guess_and_bound": midpoint_guess_and_bound,
     "rand_guess_and_bound": rand_guess_and_bound,
     "rand_distributed": rand_distributed,
     "midpoint": simple_midpoint,
 }
 
 
-def reinitialize_variables(model, config):
+def reinitialize_variables(model, config, sampler):
     """Reinitializes all variable values in the model.
 
     Excludes fixed, noncontinuous, and unbounded variables.
 
     """
-    for var in model.component_data_objects(ctype=Var, descend_into=True):
+
+    eligible_vars = []
+    for var in model._vars_list:
         if var.is_fixed() or not var.is_continuous():
             continue
         if var.lb is None or var.ub is None:
@@ -75,8 +92,34 @@ def reinitialize_variables(model, config):
                     'suppress_unbounded_warning flag.' % (var.name, var.lb, var.ub)
                 )
             continue
-        val = var.value if var.value is not None else (var.lb + var.ub) / 2
-        # apply reinitialization strategy to variable
-        var.set_value(
-            strategies[config.strategy](val, var.lb, var.ub), skip_validation=True
-        )
+
+        eligible_vars.append(var)
+
+    if config.strategy == "rand_vector":
+
+        if len(eligible_vars) == 0:
+            raise ValueError(
+                "No eligible variables to reinitialize." "Please add bounds."
+            )
+        # Collect lower and upper bounds for sampler
+        lowers = [v.lb for v in eligible_vars]
+        uppers = [v.ub for v in eligible_vars]
+
+        samples = rand_vector(lowers, uppers, sampler)
+
+        # assign samples to variables
+        for var, sample in zip(eligible_vars, samples):
+            var.set_value(sample, skip_validation=True)
+
+        return
+
+    # Otherwise
+    else:
+        for var in eligible_vars:
+            val = var.value if var.value is not None else (var.lb + var.ub) / 2
+            # print(f"val = {val}\n")
+            # apply reinitialization strategy to variable
+            var.set_value(
+                strategies[config.strategy](val, var.lb, var.ub, sampler),
+                skip_validation=True,
+            )

pyomo/contrib/multistart/tests/test_multi.py

@@ -13,6 +13,7 @@
 from io import StringIO
 
 import pyomo.common.unittest as unittest
+from pyomo.common.dependencies import attempt_import, numpy_available, scipy_available
 from pyomo.common.log import LoggingIntercept
 from pyomo.contrib.multistart.high_conf_stop import should_stop
 from pyomo.contrib.multistart.reinit import strategies
@@ -21,14 +22,18 @@
     Constraint,
     NonNegativeReals,
     Objective,
-    SolverFactory,
     Var,
     maximize,
     sin,
     value,
 )
+from pyomo.contrib.solver.common.factory import SolverFactory
+from pyomo.contrib.solver.solvers.ipopt import Ipopt
+from pyomo.contrib.solver.common.util import NoOptimalSolutionError
+from pyomo.contrib.solver.common.results import TerminationCondition
 
 
+@unittest.skipIf(not numpy_available, "Numpy not available")
 @unittest.skipIf(not SolverFactory('ipopt').available(), "IPOPT not available")
 class MultistartTests(unittest.TestCase):
     """
@@ -94,22 +99,47 @@ def test_missing_bounds(self):
                 "variable x with bounds (0, None).",
                 output.getvalue().strip(),
             )
+        with self.assertRaises(ValueError):
+            SolverFactory('multistart').solve(m, strategy="rand_vector")
 
     def test_var_value_None(self):
         m = ConcreteModel()
         m.x = Var(bounds=(0, 1))
         m.obj = Objective(expr=m.x)
         SolverFactory('multistart').solve(m)
 
+    def test_no_obj(self):
+        m = ConcreteModel()
+        m.x = Var()
+        output = StringIO()
+        with LoggingIntercept(output, 'pyomo.contrib.multistart', logging.WARNING):
+            try:
+                SolverFactory('multistart').solve(m)
+            except:
+                pass
+            self.assertIn(
+                "No objective found in the provided model. The solver will "
+                "stop if it finds a feasible solution before completing "
+                "the total number of iterations.",
+                output.getvalue().strip(),
+            )
+
     def test_model_infeasible(self):
         m = ConcreteModel()
         m.x = Var(bounds=(0, 1))
         m.c = Constraint(expr=m.x >= 2)
         m.o = Objective(expr=m.x)
-        SolverFactory('multistart').solve(m, iterations=2)
+
+        with self.assertRaises(NoOptimalSolutionError):
+            SolverFactory('multistart').solve(m, iterations=2)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.contrib.multistart', logging.WARNING):
-            SolverFactory('multistart').solve(m, iterations=-1, HCS_max_iterations=3)
+            SolverFactory('multistart').solve(
+                m,
+                iterations=-1,
+                HCS_max_iterations=3,
+                raise_exception_on_nonoptimal_result=False,
+            )
             self.assertIn(
                 "High confidence stopping rule was unable to "
                 "complete after 3 iterations.",
@@ -134,18 +164,84 @@ def test_multiple_obj(self):
         with self.assertRaisesRegex(RuntimeError, "multiple active objectives"):
             SolverFactory('multistart').solve(m)
 
-    def test_no_obj(self):
+    def test_unsupported_sampling_method(self):
         m = ConcreteModel()
-        m.x = Var()
-        with self.assertRaisesRegex(RuntimeError, "no active objective"):
-            SolverFactory('multistart').solve(m)
+        m.x = Var(bounds=(0, 1))
+        m.obj = Objective(expr=m.x)
+        with self.assertRaises(ValueError):
+            SolverFactory('multistart').solve(m, sampling_method="dummy")
 
-    def test_const_obj(self):
+    def test_unsupported_solver_type(self):
         m = ConcreteModel()
-        m.x = Var()
-        m.o = Objective(expr=5)
-        with self.assertRaisesRegex(RuntimeError, "constant objective"):
-            SolverFactory('multistart').solve(m)
+        m.x = Var(bounds=(0, 1))
+        m.obj = Objective(expr=m.x)
+        bad_solver = {"solver": "dummy"}
+        # Provided during init
+        with self.assertRaises(TypeError):
+            SolverFactory('multistart', subsolver=bad_solver).solve(m)
+        # Provided with config
+        with self.assertRaises(TypeError):
+            solver = SolverFactory('multistart')
+            solver.config.subsolver = bad_solver
+            solver.solve(m)
+        # Provided during solve
+        with self.assertRaises(TypeError):
+            solver = SolverFactory('multistart')
+            solver.config.subsolver = bad_solver
+            solver.solve(m)
+
+    def test_solver_object_matches_solver_string(self):
+        nlp_solver = Ipopt()
+        solver_str = "ipopt"
+        seed = 145
+
+        fresh_model = build_model()
+
+        m1 = fresh_model.clone()
+        results_obj_obj = SolverFactory('multistart', subsolver=nlp_solver).solve(
+            m1, seed=seed
+        )
+        m2 = fresh_model.clone()
+        results_obj_str = SolverFactory('multistart', subsolver=solver_str).solve(
+            m2, seed=seed
+        )
+        self.assertAlmostEqual(
+            results_obj_obj.incumbent_objective, results_obj_str.incumbent_objective
+        )
+
+    @unittest.skipIf(not scipy_available, "Scipy not available")
+    def test_sampling_methods(self):
+        sampling_methods = ["uniform", "latin_hypercube", "sobol"]
+        strategies = ["rand", "rand_vector"]
+        seed = 145
+        # For the simple model, check all sampling methods converge
+        simple_model = build_model()
+        for method, strategy in product(sampling_methods, strategies):
+            m = simple_model.clone()
+            res = SolverFactory("multistart").solve(
+                m, strategy=strategy, sampling_method=method, seed=seed
+            )
+            self.assertEqual(
+                res.termination_condition,
+                TerminationCondition.convergenceCriteriaSatisfied,
+            )
+
+    def test_max_time_limit(self):
+        simple_model = build_model()
+        output = StringIO()
+
+        with LoggingIntercept(output, 'pyomo.contrib.multistart', logging.WARNING):
+            res = SolverFactory('multistart').solve(
+                simple_model,
+                raise_exception_on_nonoptimal_result=False,
+                time_limit=1e-6,
+            )
+            self.assertIn(
+                "Time limit reached after 1 iterations.", output.getvalue().strip()
+            )
+            self.assertEqual(
+                res.termination_condition, TerminationCondition.maxTimeLimit
+            )
 
 
 def build_model():

pyomo/contrib/solver/tests/solvers/test_solvers.py

@@ -48,6 +48,8 @@
 from pyomo.core.expr.compare import assertExpressionsEqual
 from pyomo.core.expr.numeric_expr import LinearExpression
 
+from pyomo.contrib.multistart.multi import MultiStart
+
 np, numpy_available = attempt_import('numpy')
 parameterized, param_available = attempt_import('parameterized')
 parameterized = parameterized.parameterized
@@ -80,6 +82,7 @@ def param_as_standalone_func(cls, p, func, name):
     ('scip_persistent', ScipPersistent),
     ('gams', GAMS),
     ('knitro_direct', KnitroDirectSolver),
+    ('multistart', MultiStart),
 ]
 mip_solvers = [
     ('gurobi_persistent', GurobiPersistent),
@@ -127,6 +130,7 @@ def param_as_standalone_func(cls, p, func, name):
     ('scip_direct', ScipDirect),
     ('scip_persistent', ScipPersistent),
 ]
+meta_solvers = [('multistart', MultiStart)]
 
 
 def _load_tests(solver_list):
@@ -1165,7 +1169,7 @@ def test_results_infeasible(
         self.assertAlmostEqual(m.y.value, None)
         self.assertTrue(res.incumbent_objective is None)
 
-        if not isinstance(opt, Ipopt):
+        if not isinstance(opt, Ipopt) and not isinstance(opt, MultiStart):
             # ipopt can return the values of the variables/duals at the last iterate
             # even if it did not converge; raise_exception_on_nonoptimal_result
             # is set to False, so we are free to load infeasible solutions

pyomo/devel/initialization/__init__.py

@@ -11,4 +11,5 @@
     initialize_with_LP_approximation,
     initialize_with_piecewise_linear_approximation,
     initialize_with_global_opt,
+    initialize_with_multistart_opt,
 )

pyomo/devel/initialization/examples/init_polynomial_ex.py

@@ -53,8 +53,28 @@ def global_init_ex():
     return results.solution_status, m.x.value
 
 
+def multistart_init_ex():
+    m = build_model()
+    nlp_solver = SolverFactory('ipopt')
+    multistart_solver = SolverFactory('multistart')
+    multistart_solver.config.subsolver = nlp_solver
+
+    # multistart_solver.config.sampling_method = "uniform"
+    # multistart_solver.config.sampling_method = "lhs"
+    # multistart_solver.config.sampling_method = "sobol
+
+    multistart_solver.config.iterations = 10
+    multistart_solver.config.break_on_solution = True
+
+    results = ini.initialize_with_multistart_opt(
+        nlp=m, nlp_solver=nlp_solver, multistart_solver=multistart_solver, seed=145
+    )
+    return results, m.x.value
+
+
 if __name__ == '__main__':
     # stat, x = lp_init_ex()
     # stat, x = pwl_init_ex()
-    stat, x = global_init_ex()
+    # stat, x = global_init_ex()
+    stat, x = multistart_init_ex()
     print(stat, round(x, 4))