Adding method for multistart initialization of nonlinear models.

pyomo/contrib/multistart/high_conf_stop.py

@@ -17,6 +17,7 @@
 
 from collections import Counter
 from math import log, sqrt
+import logging
 
 
 def num_one_occurrences(observed_obj_vals, tolerance):
@@ -59,4 +60,7 @@ def should_stop(solutions, stopping_mass, stopping_delta, tolerance):
     d = stopping_delta
     c = stopping_mass
     confidence = f / n + (2 * sqrt(2) + sqrt(3)) * sqrt(log(3 / d) / n)
+    # Add temporary logger
+    logging.info(f"Number of solutions [n]:{n}; Optima viewed once [f]:{f}; \
+    Confidence:{confidence}")
     return confidence < c

pyomo/contrib/multistart/multi.py

@@ -17,11 +17,15 @@
     document_kwargs_from_configdict,
 )
 from pyomo.common.modeling import unique_component_name
+from pyomo.common.dependencies import numpy as np
 from pyomo.contrib.multistart.high_conf_stop import should_stop
 from pyomo.contrib.multistart.reinit import reinitialize_variables, strategies
 from pyomo.core import Objective, Var, minimize, value
 from pyomo.opt import SolverFactory, SolverStatus
-from pyomo.opt import TerminationCondition as tc
+from pyomo.contrib.solver.common.factory import SolverFactory as NewSolverFactory
+from pyomo.contrib.solver.common.results import SolutionStatus
+from pyomo.common.dependencies.scipy import stats
+from pyomo.common.dependencies import numpy as np
 
 logger = logging.getLogger('pyomo.contrib.multistart')
 
@@ -59,7 +63,11 @@ class MultiStart:
     )
     CONFIG.declare(
         "solver",
-        ConfigValue(default="ipopt", description="solver to use, defaults to ipopt"),
+        ConfigValue(
+            default="ipopt",
+            description="solver to use, defaults to ipopt"
+            "Should also be able to accept solver objects. In progress",
+        ),
     )
     CONFIG.declare(
         "solver_args",
@@ -120,6 +128,44 @@ class MultiStart:
             description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.",
         ),
     )
+    CONFIG.declare(
+        "break_on_solution",
+        ConfigValue(
+            default=False,
+            description="Condition to break if a feasible or optimal solution is found. Defaults to False.",
+        ),
+    )
+    CONFIG.declare(
+        "sampling_method",
+        ConfigValue(
+            default="random_uniform",
+            description="Method for sampling random starting points for reinitialization step. Supported options are \
+            'random_uniform', 'latin_hypercube', and 'sobol_sampling'",
+        ),
+    )
+    CONFIG.declare(
+        "seed",
+        ConfigValue(
+            default=None,
+            description="Seed for reproducibility in random sampling methods.",
+        ),
+    )
+    CONFIG.declare(
+        "rng",
+        ConfigValue(
+            default=None,
+            description="Random number generator for reproducibility in random sampling methods. \
+                Preferred over seed.",
+        ),
+    )
+
+    CONFIG.declare(
+        "initialize",
+        ConfigValue(
+            default=False,
+            description="Boolean for whether solver is being used to initialize model. Default is False.",
+        ),
+    )
 
     def available(self, exception_flag=True):
         """Check if solver is available.
@@ -138,8 +184,16 @@ def solve(self, model, **kwds):
         config = self.CONFIG(kwds.pop('options', {}))
         config.set_value(kwds)
 
-        # initialize the solver
-        solver = SolverFactory(config.solver)
+        # Create centralized sampler once
+        sampler = SamplingManager(
+            method=config.sampling_method, rng=config.rng, seed=config.seed
+        )
+
+        if config.initialize == True:
+            config.solver_args["load_solutions"] = False
+            config.solver_args["raise_exception_on_nonoptimal_result"] = False
+
+        solver = NewSolverFactory(config.solver)
 
         # Model sense
         objectives = model.component_data_objects(Objective, active=True)
@@ -149,14 +203,14 @@ def solve(self, model, **kwds):
             raise RuntimeError(
                 "Multistart solver is unable to handle model with multiple active objectives."
             )
-        if obj is None:
-            raise RuntimeError(
-                "Multistart solver is unable to handle model with no active objective."
-            )
-        if obj.polynomial_degree() == 0:
-            raise RuntimeError(
-                "Multistart solver received model with constant objective"
-            )
+        # if obj is None:
+        #     raise RuntimeError(
+        #         "Multistart solver is unable to handle model with no active objective."
+        #     )
+        # if obj.polynomial_degree() == 0:
+        #     raise RuntimeError(
+        #         "Multistart solver received model with constant objective"
+        #     )
 
         # store objective values and objective/result information for best
         # solution obtained
@@ -176,13 +230,22 @@ def solve(self, model, **kwds):
             )
 
             best_result = result = solver.solve(model, **config.solver_args)
-            if (
-                result.solver.status is SolverStatus.ok
-                and result.solver.termination_condition is tc.optimal
-            ):
+            # Check the solution status before loading variables into the model.
+            if config.initialize:
+                if result.solution_status in {
+                    SolutionStatus.feasible,
+                    SolutionStatus.optimal,
+                }:
+                    result.solution_loader.load_vars()
+                    logger.info(
+                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                    )
+
+            if best_result.solution_status is SolutionStatus.optimal:
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)
+
             num_iter = 0
             max_iter = config.iterations
             # if HCS rule is specified, reinitialize completely randomly until
@@ -204,15 +267,28 @@ def solve(self, model, **kwds):
                 ):
                     HCS_completed = True
                     break
+                logger.info(f"num_iter: {num_iter}\n")
                 num_iter += 1
                 # at first iteration, solve the originally passed model
                 m = model.clone() if num_iter > 1 else model
-                reinitialize_variables(m, config)
-                result = solver.solve(m, **config.solver_args)
-                if (
-                    result.solver.status is SolverStatus.ok
-                    and result.solver.termination_condition is tc.optimal
-                ):
+                reinitialize_variables(m, config, sampler)
+                result = solver.solve(m, **config.solver_args)  # , tee=True)
+
+                # Check the solution status before loading variables into the model.
+                if config.initialize:
+                    if result.solution_status in {
+                        SolutionStatus.feasible,
+                        SolutionStatus.optimal,
+                    }:
+                        result.solution_loader.load_vars()
+                        logger.info(
+                            f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                        )
+                        # If we are looking for the first feasible solution, then return immediately
+                        if config.break_on_solution:
+                            return best_result
+
+                if best_result.solution_status is SolutionStatus.optimal:
                     model_objectives = m.component_data_objects(Objective, active=True)
                     mobj = next(model_objectives)
                     obj_val = value(mobj.expr)
@@ -222,6 +298,7 @@ def solve(self, model, **kwds):
                         best_objective = obj_val
                         best_model = m
                         best_result = result
+
                 if num_iter == 1:
                     # if it's the first iteration, set the best_model and
                     # best_result regardless of solution status in case the
@@ -258,3 +335,57 @@ def __enter__(self):
 
     def __exit__(self, t, v, traceback):
         pass
+
+
+# Sampling class to organize and configure random samplers
+
+
+class SamplingManager:
+    def __init__(self, method="uniform", rng=None, seed=None):
+        aliases = {
+            "random_uniform": "uniform",
+            "uniform": "uniform",
+            "latin_hypercube": "lhs",
+            "lhs": "lhs",
+            "sobol_sampling": "sobol",
+            "sobol": "sobol",
+        }
+        self.method = aliases[method.lower()]
+
+        self.seed = seed
+
+        # Define or create a random number generator
+        # All
+
+        if rng is not None:
+            self.rng = rng
+        else:
+            self.rng = np.random.default_rng(seed)
+
+        self.qmc_sampler = None
+
+    def _ensure_qmc(self, dim):
+        if self.qmc_sampler is not None:
+            return
+
+        if self.method == "lhs":
+            self.qmc_sampler = stats.qmc.LatinHypercube(d=dim, seed=self.seed)
+        elif self.method == "sobol":
+            self.qmc_sampler = stats.qmc.Sobol(d=dim, scramble=True, seed=self.seed)
+        else:
+            raise ValueError(f"QMC sampler not valid for method '{self.method}'")
+
+    def sample_vector(self, lower, upper):
+        """Vector sample for uniform/lhs/sobol over all vars at once."""
+        lower = np.asarray(lower, dtype=float)
+        upper = np.asarray(upper, dtype=float)
+
+        if self.method == "uniform":
+            return self.rng.uniform(lower, upper)
+
+        if self.method in ("lhs", "sobol"):
+            self._ensure_qmc(dim=len(lower))
+            x = self.qmc_sampler.random(n=1)  # shape (1, d)
+            return stats.qmc.scale(x, lower, upper)[0]
+
+        raise ValueError(f"Unknown sampling method '{self.method}'")

pyomo/contrib/multistart/reinit.py

@@ -11,35 +11,39 @@
 
 import logging
 import random
+from pyomo.common.dependencies import numpy as np
+from pyomo.common.dependencies.scipy import stats
+from pyomo.core.expr.visitor import identify_variables
 
 from pyomo.core import Var
 
 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, rng):
+    sample = rng.uniform(lb, ub)  # uniform distribution between lb and ub
+    return sample
 
 
-def midpoint_guess_and_bound(val, lb, ub):
+def midpoint_guess_and_bound(val, lb, ub, rng=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, rng):
     """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)
+    return rng.uniform(val, far_bound)
 
 
-def rand_distributed(val, lb, ub, divisions=9):
+def rand_distributed(val, lb, ub, rng, 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 rng.choice(set_distributed_vals)
 
 
-def simple_midpoint(val, lb, ub):
+def simple_midpoint(val, lb, ub, rng=None):
     return (lb + ub) * 0.5
 
 
@@ -57,12 +61,15 @@ def linspace(lower, upper, n):
 }
 
 
-def reinitialize_variables(model, config):
+def reinitialize_variables(model, config, sampler):
     """Reinitializes all variable values in the model.
 
     Excludes fixed, noncontinuous, and unbounded variables.
 
     """
+
+    eligible_vars = []
+
     for var in model.component_data_objects(ctype=Var, descend_into=True):
         if var.is_fixed() or not var.is_continuous():
             continue
@@ -75,8 +82,35 @@ def reinitialize_variables(model, config):
                     'suppress_unbounded_warning flag.' % (var.name, var.lb, var.ub)
                 )
             continue
+
+        eligible_vars.append(var)
+
+    # Sample for new methods as a vector
+    if sampler.method in {"uniform", "lhs", "sobol"}:
+        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]
+
+        # Generate vector of samples using sampler
+        samples = sampler.sample_vector(lowers, uppers)
+
+        # assign samples to variables
+        for var, sample in zip(eligible_vars, samples):
+            var.set_value(sample, skip_validation=True)
+
+        return
+
+    # Otherwise use strategies to maintain original functionality
+    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), skip_validation=True
+            strategies[config.strategy](val, var.lb, var.ub, sampler.rng),
+            skip_validation=True,
         )

pyomo/contrib/multistart/tests/test_multi.py

@@ -134,18 +134,19 @@ def test_multiple_obj(self):
         with self.assertRaisesRegex(RuntimeError, "multiple active objectives"):
             SolverFactory('multistart').solve(m)
 
-    def test_no_obj(self):
-        m = ConcreteModel()
-        m.x = Var()
-        with self.assertRaisesRegex(RuntimeError, "no active objective"):
-            SolverFactory('multistart').solve(m)
-
-    def test_const_obj(self):
-        m = ConcreteModel()
-        m.x = Var()
-        m.o = Objective(expr=5)
-        with self.assertRaisesRegex(RuntimeError, "constant objective"):
-            SolverFactory('multistart').solve(m)
+    # Would like to remove these tests to allow for square model solves.
+    # def test_no_obj(self):
+    #     m = ConcreteModel()
+    #     m.x = Var()
+    #     with self.assertRaisesRegex(RuntimeError, "no active objective"):
+    #         SolverFactory('multistart').solve(m)
+
+    # def test_const_obj(self):
+    #     m = ConcreteModel()
+    #     m.x = Var()
+    #     m.o = Objective(expr=5)
+    #     with self.assertRaisesRegex(RuntimeError, "constant objective"):
+    #         SolverFactory('multistart').solve(m)
 
 
 def build_model():