diff --git a/pyomo/contrib/multistart/high_conf_stop.py b/pyomo/contrib/multistart/high_conf_stop.py index b9efdfd65f9..91fd0244775 100644 --- a/pyomo/contrib/multistart/high_conf_stop.py +++ b/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 diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py index 32fb8ce3cff..92e0f95d8cc 100644 --- a/pyomo/contrib/multistart/multi.py +++ b/pyomo/contrib/multistart/multi.py @@ -12,23 +12,208 @@ from pyomo.common.config import ( ConfigBlock, + ConfigDict, ConfigValue, In, document_kwargs_from_configdict, + document_class_CONFIG, + document_configdict, + ADVANCED_OPTION, ) + +from typing import Any, Optional +import datetime +from pyomo.common.timing import HierarchicalTimer, default_timer 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.core import Objective, Constraint, Var, minimize, value +from pyomo.contrib.solver.common.base import SolverBase +from pyomo.contrib.solver.common.config import SolverConfig +from pyomo.contrib.solver.common.factory import SolverFactory +from pyomo.contrib.solver.common.results import ( + Results, + TerminationCondition, + SolutionStatus, +) + +from pyomo.contrib.solver.common.util import ( + NoOptimalSolutionError, + NoSolutionError, + NoFeasibleSolutionError, +) +from pyomo.util.vars_from_expressions import get_vars_from_components + +from pyomo.common.dependencies.scipy import stats +from pyomo.common.dependencies import numpy as np +from pyomo.core.staleflag import StaleFlagManager logger = logging.getLogger('pyomo.contrib.multistart') +@document_configdict() +class MultistartConfig(SolverConfig): + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + super().__init__( + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + + self.strategy = self.declare( + "strategy", + ConfigValue( + default="rand", + domain=In(strategies.keys()), + description="Specify the restart strategy. Defaults to rand.", + doc="""Specify the restart strategy. + + - "rand": random choice between variable bounds + - "rand_vector": random choice, vectorized approach with sampler + - "midpoint_guess_and_bound": midpoint between current value and farthest bound + - "rand_guess_and_bound": random choice between current value and farthest bound + - "rand_distributed": random choice among evenly distributed values + - "midpoint": exact midpoint between the bounds. If using this option, multiple iterations are useless. + """, + ), + ) + self.subsolver = self.declare( + "subsolver", + ConfigValue( + default="ipopt", + description="solver to use, defaults to ipopt" + "Should also be able to accept solver objects. In progress", + ), + ) + self.subsolver_args = self.declare( + "subsolver_args", + ConfigValue( + default={}, + description="Dictionary of keyword arguments to pass to the solver.", + ), + ) + self.iterations = self.declare( + "iterations", + ConfigValue( + default=10, + description="Specify the number of iterations, defaults to 10. " + "If -1 is specified, the high confidence stopping rule will be used", + ), + ) + self.stopping_mass = self.declare( + "stopping_mass", + ConfigValue( + default=0.5, + description="Maximum allowable estimated missing mass of optima.", + doc="""Maximum allowable estimated missing mass of optima for the + high confidence stopping rule, only used with the random strategy. + The lower the parameter, the stricter the rule. + Value bounded in (0, 1].""", + ), + ) + self.stopping_delta = self.declare( + "stopping_delta", + ConfigValue( + default=0.5, + description="1 minus the confidence level required for the stopping rule.", + doc="""1 minus the confidence level required for the stopping rule for the + high confidence stopping rule, only used with the random strategy. + The lower the parameter, the stricter the rule. visibility=DEVELOPER_OPTION, + Value bounded in (0, 1].""", + ), + ) + self.suppress_unbounded_warning = self.declare( + "suppress_unbounded_warning", + ConfigValue( + default=False, + domain=bool, + description="True to suppress warning for skipping unbounded variables.", + ), + ) + self.HCS_max_iterations = self.declare( + "HCS_max_iterations", + ConfigValue( + default=1000, + description="Maximum number of iterations before interrupting the high confidence stopping rule.", + ), + ) + self.HCS_tolerance = self.declare( + "HCS_tolerance", + ConfigValue( + default=0, + description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.", + ), + ) + self.break_on_solution = self.declare( + "break_on_solution", + ConfigValue( + default=False, + description="Condition to break if a feasible or optimal solution is found. Defaults to False.", + ), + ) + self.sampling_method = self.declare( + "sampling_method", + ConfigValue( + default="latin_hypercube", + description="Method for sampling random starting points for reinitialization step. " + "Supported options are 'random_uniform', 'latin_hypercube', and 'sobol_sampling'. " + "Only utilized when config.strategy is 'rand_vector'.", + ), + ) + self.seed = self.declare( + "seed", + ConfigValue( + default=None, + description="Seed for reproducibility in random sampling methods.", + ), + ) + self.rng = self.declare( + "rng", + ConfigValue( + default=None, + description="Random number generator for reproducibility in random sampling methods. \ + Preferred over seed.", + ), + ) + + +class MultiStartResults(Results): + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + super().__init__( + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + self.feasible_solution_list: Optional[list] = self.declare( + 'feasible_solution_list', + ConfigValue( + description="Object for loading the solution back into the model." + ), + ) + + @SolverFactory.register('multistart', doc='MultiStart solver for NLPs') -@document_kwargs_from_configdict('CONFIG') -class MultiStart: +@document_class_CONFIG(methods=['solve']) +class MultiStart(SolverBase): """Solver wrapper that initializes at multiple starting points. # TODO: also return appropriate duals @@ -40,221 +225,284 @@ class MultiStart: """ - CONFIG = ConfigBlock("MultiStart") - CONFIG.declare( - "strategy", - ConfigValue( - default="rand", - domain=In(strategies.keys()), - description="Specify the restart strategy. Defaults to rand.", - doc="""Specify the restart strategy. - - - "rand": random choice between variable bounds - - "midpoint_guess_and_bound": midpoint between current value and farthest bound - - "rand_guess_and_bound": random choice between current value and farthest bound - - "rand_distributed": random choice among evenly distributed values - - "midpoint": exact midpoint between the bounds. If using this option, multiple iterations are useless. - """, - ), - ) - CONFIG.declare( - "solver", - ConfigValue(default="ipopt", description="solver to use, defaults to ipopt"), - ) - CONFIG.declare( - "solver_args", - ConfigValue( - default={}, - description="Dictionary of keyword arguments to pass to the solver.", - ), - ) - CONFIG.declare( - "iterations", - ConfigValue( - default=10, - description="Specify the number of iterations, defaults to 10. " - "If -1 is specified, the high confidence stopping rule will be used", - ), - ) - CONFIG.declare( - "stopping_mass", - ConfigValue( - default=0.5, - description="Maximum allowable estimated missing mass of optima.", - doc="""Maximum allowable estimated missing mass of optima for the - high confidence stopping rule, only used with the random strategy. - The lower the parameter, the stricter the rule. - Value bounded in (0, 1].""", - ), - ) - CONFIG.declare( - "stopping_delta", - ConfigValue( - default=0.5, - description="1 minus the confidence level required for the stopping rule.", - doc="""1 minus the confidence level required for the stopping rule for the - high confidence stopping rule, only used with the random strategy. - The lower the parameter, the stricter the rule. - Value bounded in (0, 1].""", - ), - ) - CONFIG.declare( - "suppress_unbounded_warning", - ConfigValue( - default=False, - domain=bool, - description="True to suppress warning for skipping unbounded variables.", - ), - ) - CONFIG.declare( - "HCS_max_iterations", - ConfigValue( - default=1000, - description="Maximum number of iterations before interrupting the high confidence stopping rule.", - ), - ) - CONFIG.declare( - "HCS_tolerance", - ConfigValue( - default=0, - description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.", - ), - ) + CONFIG = MultistartConfig() + + def __init__(self, **kwds: Any) -> None: + super().__init__(**kwds) + + #: Instance configuration; + self.config = self.config def available(self, exception_flag=True): """Check if solver is available. - TODO: For now, it is always available. However, sub-solvers may not - always be available, and so this should reflect that possibility. + The multistart solver wrapper should always be available, + but it is not guaranteed the subsolvers will be. + Check if the selected subsolver is available, which by default is ipopt. + """ + + subsolver = SolverFactory(self.config.subsolver) + return subsolver.available() + def version(self): + """Get solver version.""" """ - return True + Original implementation: 0.1.0, + Current implementation: 0.2.0, + """ + current = (0, 2, 0) + return current def license_is_valid(self): return True def solve(self, model, **kwds): + start_time = default_timer() + # initialize keyword args - config = self.CONFIG(kwds.pop('options', {})) + config = self.config(kwds.pop('options', {})) config.set_value(kwds) - # initialize the solver - solver = SolverFactory(config.solver) + timer = config.timer + if timer is None: + timer = config.timer = HierarchicalTimer() + + # Allocate the results object so we can populate it as we go + results = MultiStartResults() + results.timing_info.start_timestamp = datetime.datetime.now( + datetime.timezone.utc + ) + # As we are about to run a solver, update the stale flag + StaleFlagManager.mark_all_as_stale() + + # Create centralized sampler once + sampler = SamplingManager( + method=config.sampling_method, rng=config.rng, seed=config.seed + ) + + # Set sub-solver options + config.subsolver_args["load_solutions"] = False + config.subsolver_args["raise_exception_on_nonoptimal_result"] = False + config.subsolver_args["time_limit"] = config.time_limit + config.subsolver_args["tee"] = config.tee + + solver = SolverFactory(config.subsolver) # Model sense - objectives = model.component_data_objects(Objective, active=True) - obj = next(objectives, None) - # Check model validity - if next(objectives, None) is not None: + objectives = list(model.component_data_objects(Objective, active=True)) + # Check length + if len(objectives) > 1: 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" - ) + elif len(objectives) == 1: + obj = objectives[0] + obj.sign = 1 if obj.sense == minimize else -1 + obj_sign = obj.sign + + else: + obj = None + obj_sign = 1 + config.break_on_solution = True # store objective values and objective/result information for best # solution obtained objectives = [] - obj_sign = 1 if obj.sense == minimize else -1 - best_objective = float('inf') * obj_sign - best_model = model best_result = None + best_objective = float('inf') * obj_sign + results.feasible_solution_list = [] - try: - # create temporary variable list for value transfer - tmp_var_list_name = unique_component_name(model, "_vars_list") + timer.start('initial_solve') + # create temporary variable list for value transfer + tmp_var_list_name = unique_component_name(model, "_vars_list") + setattr( + model, + tmp_var_list_name, + list(model.component_data_objects(Var, descend_into=True)), + ) + # If the list has nothing in it, check components + if len(model._vars_list) == 0: setattr( model, tmp_var_list_name, - list(model.component_data_objects(ctype=Var, descend_into=True)), + list( + get_vars_from_components( + model, ctype=(Constraint, Objective), active=True + ) + ), + ) + best_result = result = solver.solve(model, **config.subsolver_args) + # Check the solution status before loading variables into the model. + if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: + results.feasible_solution_list.append(result) + logger.info( + f'solved NLP: {result.solution_status}, {result.termination_condition}' ) - best_result = result = solver.solve(model, **config.solver_args) - if ( - result.solver.status is SolverStatus.ok - and result.solver.termination_condition is tc.optimal - ): - obj_val = value(obj.expr) + if result.solution_status is SolutionStatus.optimal: + if obj is not None: + obj_val = result.incumbent_objective best_objective = obj_val objectives.append(obj_val) - num_iter = 0 - max_iter = config.iterations - # if HCS rule is specified, reinitialize completely randomly until - # rule specifies stopping - using_HCS = config.iterations == -1 - HCS_completed = False - if using_HCS: - assert ( - config.strategy == "rand" - ), "High confidence stopping rule requires rand strategy." - max_iter = config.HCS_max_iterations - - while num_iter < max_iter: - if using_HCS and should_stop( - objectives, - config.stopping_mass, - config.stopping_delta, - config.HCS_tolerance, - ): - HCS_completed = True + timer.stop('initial_solve') + + num_iter = 0 + max_iter = config.iterations + # if HCS rule is specified, reinitialize completely randomly until + # rule specifies stopping + using_HCS = config.iterations == -1 + HCS_completed = False + if using_HCS: + assert ( + config.strategy == "rand" + ), "High confidence stopping rule requires rand strategy." + max_iter = config.HCS_max_iterations + + timer.start('iterative_solves') + while num_iter < max_iter: + # timer.start(f"timer_iter_{num_iter}") + if using_HCS and should_stop( + objectives, + config.stopping_mass, + config.stopping_delta, + config.HCS_tolerance, + ): + HCS_completed = True + # timer.stop(f"timer_iter_{num_iter}") + break + + num_iter += 1 + logger.info(f"num_iter: {num_iter}\n") + + # at first iteration, solve the originally passed model + m = model + reinitialize_variables(m, config, sampler) + result = solver.solve(m, **config.subsolver_args) + # if config.load_solutions: + # Check the solution status before loading variables into the model. + if result.solution_status in { + SolutionStatus.feasible, + SolutionStatus.optimal, + }: + logger.info( + f'solved NLP: {result.solution_status}, {result.termination_condition}' + ) + results.feasible_solution_list.append(result) + # If we are looking for the first feasible solution, then return immediately + if config.break_on_solution: + best_result = result + # timer.stop(f"timer_iter_{num_iter}") break - 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 - ): - model_objectives = m.component_data_objects(Objective, active=True) - mobj = next(model_objectives) - obj_val = value(mobj.expr) + + if result.solution_status is SolutionStatus.optimal: + if obj is not None: + obj_val = result.incumbent_objective objectives.append(obj_val) if obj_val * obj_sign < obj_sign * best_objective: # objective has improved 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 - # model is infeasible. - best_model = m - best_result = result + # timer.stop(f"timer_iter_{num_iter}") - if using_HCS and not HCS_completed: + timer.stop('iterative_solves') + delattr(model, tmp_var_list_name) + + if using_HCS: + if not HCS_completed: logger.warning( "High confidence stopping rule was unable to complete " "after %s iterations. To increase this limit, change the " "HCS_max_iterations flag." % num_iter ) - # if no better result was found than initial solve, then return - # that without needing to copy variables. - if best_model is model: - return best_result + if ( + config.raise_exception_on_nonoptimal_result + and best_result.solution_status != SolutionStatus.optimal + ): + raise NoOptimalSolutionError() + + results.solution_loader = best_result.solution_loader + results.termination_condition = best_result.termination_condition + results.solution_status = best_result.solution_status + results.incumbent_objective = best_result.incumbent_objective + results.solver_log = best_result.solver_log + + if config.load_solutions: + if results.solution_status == SolutionStatus.noSolution: + raise NoSolutionError() - # reassign the given models vars to the new models vars - orig_var_list = getattr(model, tmp_var_list_name) - best_soln_var_list = getattr(best_model, tmp_var_list_name) - for orig_var, new_var in zip(orig_var_list, best_soln_var_list): - if not orig_var.is_fixed(): - orig_var.set_value(new_var.value, skip_validation=True) + results.solution_loader.load_solution() - return best_result - finally: - # Remove temporary variable list - delattr(model, tmp_var_list_name) + results.solver_name = self.name + results.solver_version = self.version() + results.solver_config = config + results.timing_info.timer = timer + results.timing_info.wall_time = default_timer() - start_time + return results def __enter__(self): return self def __exit__(self, t, v, traceback): pass + + +# Sampling class to organize and configure random samplers +class SamplingManager: + def __init__(self, method="lhs", 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 + 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, rng=self.rng) + elif self.method == "sobol": + self.qmc_sampler = stats.qmc.Sobol(d=dim, scramble=True, seed=self.rng) + else: + raise ValueError(f"QMC sampler not valid for method '{self.method}'") + + def sample_scalar(self, lower, upper): + if self.method == "uniform": + return self.rng.uniform(lower, upper) + + if self.method in ("lhs", "sobol"): + self._ensure_qmc(dim=1) + x = self.qmc_sampler.random(n=1) # shape (1, d) + return stats.qmc.scale(x, lower, upper).item() + + raise ValueError(f"Unknown sampling 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}'") diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py index a3b52a8d611..04aa40cf4a1 100644 --- a/pyomo/contrib/multistart/reinit.py +++ b/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,6 +64,7 @@ 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, @@ -57,13 +72,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. """ - 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, + ) diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py index 1c34138fbb3..cd6ca6878c2 100644 --- a/pyomo/contrib/multistart/tests/test_multi.py +++ b/pyomo/contrib/multistart/tests/test_multi.py @@ -21,12 +21,13 @@ Constraint, NonNegativeReals, Objective, - SolverFactory, Var, maximize, sin, value, ) +from pyomo.contrib.solver.common.factory import SolverFactory +from pyomo.contrib.solver.common.util import NoOptimalSolutionError @unittest.skipIf(not SolverFactory('ipopt').available(), "IPOPT not available") @@ -106,10 +107,17 @@ def test_model_infeasible(self): 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,19 +142,6 @@ 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) - def build_model(): """Simple non-convex model with many local minima""" diff --git a/pyomo/contrib/solver/common/solution_loader.py b/pyomo/contrib/solver/common/solution_loader.py index faaf7c75685..2fda01e5c32 100644 --- a/pyomo/contrib/solver/common/solution_loader.py +++ b/pyomo/contrib/solver/common/solution_loader.py @@ -63,7 +63,7 @@ def solution(self, solution_id: Any) -> "SolutionLoaderView": results = solver.solve(model) results.solution(2).load_vars() - results.solution(2).load_import_suffixes() + results.solution(l2).load_import_suffixes() Parameters ---------- diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 0b30a6eb923..9fc95e32d91 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/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), @@ -1150,7 +1153,7 @@ def test_results_infeasible( opt.config.load_solutions = False res = opt.solve(m) self.assertNotEqual(res.solution_status, SolutionStatus.optimal) - if isinstance(opt, Ipopt): + if isinstance(opt, Ipopt) or isinstance(opt, MultiStart): acceptable_termination_conditions = { TerminationCondition.locallyInfeasible, TerminationCondition.unbounded, @@ -1165,7 +1168,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 @@ -1229,7 +1232,7 @@ def test_trivial_constraints( opt.config.load_solutions = False res = opt.solve(m) self.assertNotEqual(res.solution_status, SolutionStatus.optimal) - if isinstance(opt, Ipopt): + if isinstance(opt, Ipopt) or isinstance(opt, MultiStart): acceptable_termination_conditions = { TerminationCondition.locallyInfeasible, TerminationCondition.unbounded, @@ -1909,7 +1912,7 @@ def test_time_limit( constant=0, ) m.c2[t] = expr == 1 - if isinstance(opt, Ipopt): + if isinstance(opt, Ipopt) or isinstance(opt, MultiStart): opt.config.time_limit = 1e-6 else: opt.config.time_limit = 0 diff --git a/pyomo/devel/initialization/__init__.py b/pyomo/devel/initialization/__init__.py index 49da58c9d5e..d3796e90c3f 100644 --- a/pyomo/devel/initialization/__init__.py +++ b/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, ) diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py index 643c2faf180..7402ab55321 100644 --- a/pyomo/devel/initialization/examples/init_polynomial_ex.py +++ b/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.strategy = "rand" + # multistart_solver.config.strategy = "rand_vector" + # 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)) diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py index 83b8afba98e..3f3ecb9fee0 100644 --- a/pyomo/devel/initialization/global_init.py +++ b/pyomo/devel/initialization/global_init.py @@ -31,6 +31,14 @@ def _initialize_with_global_solver( 'interfaces, so the global solvers are limited to ScipDirect, ' 'ScipPersistent, and GurobiDirectMINLP.' ) + # Check if time limit is provided for global solver + if global_solver.config.time_limit is None: + logger.warning( + 'No time limit set for global optimizer. ' + 'For a large model, this may take a long time. ' + 'Consider setting a time limit using global_solver.config.time_limit.' + ) + res = global_solver.solve( nlp, load_solutions=True, @@ -40,15 +48,5 @@ def _initialize_with_global_solver( logger.info( f'solved NLP with {global_solver.name}: {res.solution_status}, {res.termination_condition}' ) - res = nlp_solver.solve( - nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False - ) - logger.info( - f'solved NLP with {nlp_solver.name}: {res.solution_status}, {res.termination_condition}' - ) - if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: - res.solution_loader.load_vars() - else: - logger.warning('initialization was not successful via global optimization') return res diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py index 7b1117cf5b9..55862e79402 100644 --- a/pyomo/devel/initialization/initialize.py +++ b/pyomo/devel/initialization/initialize.py @@ -18,10 +18,14 @@ from pyomo.devel.initialization.lp_approx_init import _initialize_with_LP_approximation from pyomo.contrib.solver.common.base import SolverBase from pyomo.devel.initialization.global_init import _initialize_with_global_solver +from pyomo.devel.initialization.multistart_init import ( + _initialize_with_multistart_solver, +) from pyomo.contrib.solver.common.factory import SolverFactory from pyomo.contrib.solver.common.results import Results import logging from pyomo.contrib.solver.common.results import SolutionStatus +import pyomo.environ as pyo logger = logging.getLogger(__name__) @@ -77,6 +81,21 @@ def _try_nlp_solve(nlp: BlockData, nlp_solver: SolverBase): return res +def _retry_nlp_solve(nlp: BlockData, nlp_solver: SolverBase): + # retry to solve the original nlp after using an initialization method + nlp_res = nlp_solver.solve( + nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + logger.info(f'resolved NLP with {nlp_solver.name}: {nlp_res.solution_status}, \ + {nlp_res.termination_condition}') + if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: + nlp_res.solution_loader.load_vars() + else: + logger.warning('initialization did not find feasible solution') + + return nlp_res + + def initialize_with_piecewise_linear_approximation( nlp: BlockData, nlp_solver: SolverBase | None = None, @@ -156,7 +175,9 @@ def initialize_with_piecewise_linear_approximation( finally: _cleanup(orig_var_data) - return res + nlp_res = _retry_nlp_solve(nlp, nlp_solver) + + return nlp_res def initialize_with_LP_approximation( @@ -239,7 +260,9 @@ def initialize_with_LP_approximation( finally: _cleanup(orig_var_data) - return res + nlp_res = _retry_nlp_solve(nlp, nlp_solver) + + return nlp_res def initialize_with_global_opt( @@ -293,4 +316,72 @@ def initialize_with_global_opt( finally: _cleanup(orig_var_data) - return res + nlp_res = _retry_nlp_solve(nlp, nlp_solver) + + return nlp_res + + +def initialize_with_multistart_opt( + nlp: BlockData, + nlp_solver: SolverBase | None = None, + multistart_solver=None, + skip_initial_nlp_solve: bool = False, + default_bound: float = 1e8, + seed=0, +) -> Results: + """ + Attempt to initialize and subsequently solve the model given by ``nlp``. + The basic idea is to apply some method to find good initial values for + the variables and then try to solve the problem with ``nlp_solver``. + + Parameters + ---------- + nlp: BlockData + The pyomo model to be initialized. + nlp_solver: Optional[SolverBase] + A solver interface appropriate for NLPs. + Default: ipopt + multistart_solver: Optional[SolverFactory] + A configured multistart solver object for performing multistart optimization + skip_initial_nlp_solve: bool + If True, the initial attempt at solving the NLP without initialization + will be skipped. + seed: 0 + Set reproducibility seed to make result deterministic. + + Returns + ------- + res: pyomo.contrib.solver.common.results.Results + The results object obtained the last time the nlp_solver was used to + try and solve the model. + """ + + if nlp_solver is None: + nlp_solver = _get_solver('ipopt', 'local NLP solver') + + if multistart_solver is None: + multistart_solver = pyo.SolverFactory("multistart") + multistart_solver.config.seed = seed + multistart_solver.config.sampling_method = "lhs" + multistart_solver.config.break_on_solution = True + + if not skip_initial_nlp_solve: + res = _try_nlp_solve(nlp, nlp_solver) + if res.solution_status == SolutionStatus.optimal: + return res + + orig_var_data = _setup(nlp) + + try: + res = _initialize_with_multistart_solver( + nlp=nlp, + multistart_solver=multistart_solver, + seed=seed, + default_bound=default_bound, + ) + finally: + _cleanup(orig_var_data) + + nlp_res = _retry_nlp_solve(nlp, nlp_solver) + + return nlp_res diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py index c8f7c2ff4fc..d20c92c0203 100644 --- a/pyomo/devel/initialization/lp_approx_init.py +++ b/pyomo/devel/initialization/lp_approx_init.py @@ -218,17 +218,4 @@ def _initialize_with_LP_approximation( ) logger.info(f'solved LP: {lp_res.solution_status}, {lp_res.termination_condition}') - # try solving the NLP - nlp_res = nlp_solver.solve( - orig_nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False - ) - logger.info( - f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}' - ) - - if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: - nlp_res.solution_loader.load_vars() - else: - logger.warning('initialization was not successful via LP approximation') - - return nlp_res + return lp_res diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py new file mode 100644 index 00000000000..0f33c3d7e45 --- /dev/null +++ b/pyomo/devel/initialization/multistart_init.py @@ -0,0 +1,47 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from pyomo.core.base.block import BlockData +from pyomo.contrib.solver.common.base import SolverBase +from pyomo.common.dependencies import numpy as np +import pyomo.environ as pyo +from pyomo.contrib.solver.common.results import SolutionStatus +from pyomo.devel.initialization.bounds.bound_variables import ( + bound_all_nonlinear_variables, +) +from pyomo.devel.initialization.utils import shallow_clone, fix_vars_with_equal_bounds +import logging + +logger = logging.getLogger(__name__) + + +def _initialize_with_multistart_solver( + nlp: BlockData, multistart_solver, seed, default_bound=1.0e8 +): + + # Make a shallow clone + nlp = shallow_clone(nlp) + # bounds on the nonlinear variables + bound_all_nonlinear_variables(nlp, default_bound=default_bound) + logger.info('bounded nonlinear variables') + + # fix variables with equal bounds for sampler + fix_vars_with_equal_bounds(nlp) + logger.info('fixed variables with equal bounds') + + multistart_solver.config.seed = seed + multistart_solver.config.load_solutions = False + multistart_solver.config.raise_exception_on_nonoptimal_result = False + + res = multistart_solver.solve(nlp) + logger.info('Finished multistart optimization iterations.') + if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: + res.solution_loader.load_vars() + + return res diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py index 67a260d9904..1d39a9c3dae 100644 --- a/pyomo/devel/initialization/tests/test_initialization.py +++ b/pyomo/devel/initialization/tests/test_initialization.py @@ -16,6 +16,7 @@ lp_init_ex, pwl_init_ex, global_init_ex, + multistart_init_ex, ) from pyomo.common import unittest from pyomo.common.dependencies import scipy_available @@ -85,6 +86,11 @@ def test_poly_lp(self): self.assertEqual(stat, SolutionStatus.optimal) self.assertAlmostEqual(x, -9.920159607881597) + def test_poly_multistart(self): + stat, x = lp_init_ex() + self.assertEqual(stat, SolutionStatus.optimal) + self.assertAlmostEqual(x, -9.920159607881597) + class TestInit(unittest.TestCase): @unittest.skipUnless(highs.available(), 'highs is not available') @@ -196,6 +202,8 @@ def test_pwl_init(self): 25: ([-9.91992877683681], 1e-6, 1e-6), 26: ([-9.920038488200985], 1e-6, 1e-6), 27: ([-9.920096055464825], 1e-6, 1e-6), + # For new second solve, repeat last value. + 28: ([-9.920096055464825], 1e-6, 1e-6), }, ) mip_solver = SolverFactory('highs') @@ -227,8 +235,8 @@ def test_pwl_ineq(self): if __name__ == '__main__': + import logging logging.basicConfig(level=logging.INFO) t = TestInit() - t.test_pwl_init()