From a221740ac54f00ca813ff82b60f1e072d251cffe Mon Sep 17 00:00:00 2001 From: Alinson S Xavier Date: Fri, 25 Sep 2020 06:02:07 -0500 Subject: [PATCH] Implement lazy callbacks & two-phase gap --- miplearn/components/component.py | 3 + miplearn/components/lazy_static.py | 52 +++++++-- miplearn/components/tests/test_lazy_static.py | 6 +- miplearn/log.py | 47 ++++---- miplearn/solvers/gurobi.py | 108 ++++++++++++++++-- miplearn/solvers/internal.py | 19 ++- miplearn/solvers/learning.py | 86 ++++++++++---- miplearn/solvers/pyomo/base.py | 11 +- miplearn/solvers/tests/__init__.py | 4 +- miplearn/solvers/tests/test_lazy_cb.py | 27 +++++ 10 files changed, 288 insertions(+), 75 deletions(-) create mode 100644 miplearn/solvers/tests/test_lazy_cb.py diff --git a/miplearn/components/component.py b/miplearn/components/component.py index d791b4e..5c353e0 100644 --- a/miplearn/components/component.py +++ b/miplearn/components/component.py @@ -24,3 +24,6 @@ class Component(ABC): def after_iteration(self, solver, instance, model): return False + + def on_lazy_callback(self, solver, instance, model): + return diff --git a/miplearn/components/lazy_static.py b/miplearn/components/lazy_static.py index 289eb2d..d223eaf 100644 --- a/miplearn/components/lazy_static.py +++ b/miplearn/components/lazy_static.py @@ -21,14 +21,29 @@ class LazyConstraint: class StaticLazyConstraintsComponent(Component): def __init__(self, classifier=CountingClassifier(), - threshold=0.05): + threshold=0.05, + use_two_phase_gap=True, + large_gap=1e-2, + violation_tolerance=-0.5, + ): self.threshold = threshold self.classifier_prototype = classifier self.classifiers = {} self.pool = [] + self.original_gap = None + self.large_gap = large_gap + self.is_gap_large = False + self.use_two_phase_gap = use_two_phase_gap + self.violation_tolerance = violation_tolerance def before_solve(self, solver, instance, model): self.pool = [] + if not solver.use_lazy_cb and self.use_two_phase_gap: + logger.info("Increasing gap tolerance to %f", self.large_gap) + self.original_gap = solver.gap_tolerance + self.is_gap_large = True + solver.internal_solver.set_gap_tolerance(self.large_gap) + instance.found_violated_lazy_constraints = [] if instance.has_static_lazy_constraints(): self._extract_and_predict_static(solver, instance) @@ -37,23 +52,41 @@ class StaticLazyConstraintsComponent(Component): pass def after_iteration(self, solver, instance, model): - logger.info("Finding violated lazy constraints...") + if solver.use_lazy_cb: + return False + else: + should_repeat = self._check_and_add(instance, solver) + if should_repeat: + return True + else: + if self.is_gap_large: + logger.info("Restoring gap tolerance to %f", self.original_gap) + solver.internal_solver.set_gap_tolerance(self.original_gap) + self.is_gap_large = False + return True + else: + return False + + def on_lazy_callback(self, solver, instance, model): + self._check_and_add(instance, solver) + + def _check_and_add(self, instance, solver): + logger.debug("Finding violated lazy constraints...") constraints_to_add = [] for c in self.pool: - if not solver.internal_solver.is_constraint_satisfied(c.obj): + if not solver.internal_solver.is_constraint_satisfied(c.obj, + tol=self.violation_tolerance): constraints_to_add.append(c) for c in constraints_to_add: self.pool.remove(c) solver.internal_solver.add_constraint(c.obj) instance.found_violated_lazy_constraints += [c.cid] if len(constraints_to_add) > 0: - logger.info("Added %d lazy constraints back into the model" % len(constraints_to_add)) - logger.info("Lazy constraint pool has %d constraints" % len(self.pool)) + logger.info("%8d lazy constraints added %8d in the pool" % (len(constraints_to_add), len(self.pool))) return True else: - logger.info("Found no violated lazy constraints") return False - + def fit(self, training_instances): training_instances = [t for t in training_instances @@ -92,7 +125,7 @@ class StaticLazyConstraintsComponent(Component): obj=solver.internal_solver.extract_constraint(cid)) constraints[category] += [c] self.pool.append(c) - logger.info("Extracted %d lazy constraints" % len(self.pool)) + logger.info("%8d lazy constraints extracted" % len(self.pool)) logger.info("Predicting required lazy constraints...") n_added = 0 for (category, x_values) in x.items(): @@ -108,8 +141,7 @@ class StaticLazyConstraintsComponent(Component): self.pool.remove(c) solver.internal_solver.add_constraint(c.obj) instance.found_violated_lazy_constraints += [c.cid] - logger.info("Added %d lazy constraints back into the model" % n_added) - logger.info("Lazy constraint pool has %d constraints" % len(self.pool)) + logger.info("%8d lazy constraints added %8d in the pool" % (n_added, len(self.pool))) def _collect_constraints(self, train_instances): constraints = {} diff --git a/miplearn/components/tests/test_lazy_static.py b/miplearn/components/tests/test_lazy_static.py index c4b651d..da7bfb9 100644 --- a/miplearn/components/tests/test_lazy_static.py +++ b/miplearn/components/tests/test_lazy_static.py @@ -13,6 +13,9 @@ from miplearn.classifiers import Classifier def test_usage_with_solver(): solver = Mock(spec=LearningSolver) + solver.use_lazy_cb = False + solver.gap_tolerance = 1e-4 + internal = solver.internal_solver = Mock(spec=InternalSolver) internal.get_constraint_ids = Mock(return_value=["c1", "c2", "c3", "c4"]) internal.extract_constraint = Mock(side_effect=lambda cid: "<%s>" % cid) @@ -37,7 +40,8 @@ def test_usage_with_solver(): "c4": "type-b", }[cid]) - component = StaticLazyConstraintsComponent(threshold=0.90) + component = StaticLazyConstraintsComponent(threshold=0.90, + use_two_phase_gap=False) component.classifiers = { "type-a": Mock(spec=Classifier), "type-b": Mock(spec=Classifier), diff --git a/miplearn/log.py b/miplearn/log.py index c41f0cc..96f56d8 100644 --- a/miplearn/log.py +++ b/miplearn/log.py @@ -7,40 +7,43 @@ import logging import time import sys -if sys.stdout.isatty(): - log_colors = { - "green": '\033[92m', - "yellow": '\033[93m', - "red": '\033[91m', - "reset": '\033[0m', - } -else: - log_colors = { - "green": "", - "yellow": "", - "red": "", - "reset": "", - } - class TimeFormatter(): - def __init__(self, start_time): + def __init__(self, start_time, log_colors): self.start_time = start_time + self.log_colors = log_colors def format(self, record): if record.levelno >= logging.ERROR: - color = log_colors["red"] + color = self.log_colors["red"] elif record.levelno >= logging.WARNING: - color = log_colors["yellow"] + color = self.log_colors["yellow"] else: - color = log_colors["green"] + color = self.log_colors["green"] return "%s[%12.3f]%s %s" % (color, record.created - self.start_time, - log_colors["reset"], + self.log_colors["reset"], record.getMessage()) -def setup_logger(start_time): +def setup_logger(start_time=None, + force_color=False): + if start_time is None: + start_time = time.time() + if sys.stdout.isatty() or force_color: + log_colors = { + "green": '\033[92m', + "yellow": '\033[93m', + "red": '\033[91m', + "reset": '\033[0m', + } + else: + log_colors = { + "green": "", + "yellow": "", + "red": "", + "reset": "", + } handler = logging.StreamHandler() - handler.setFormatter(TimeFormatter(start_time)) + handler.setFormatter(TimeFormatter(start_time, log_colors)) logging.getLogger().addHandler(handler) logging.getLogger("miplearn").setLevel(logging.INFO) lg = logging.getLogger("miplearn") diff --git a/miplearn/solvers/gurobi.py b/miplearn/solvers/gurobi.py index 0587e0f..cae087b 100644 --- a/miplearn/solvers/gurobi.py +++ b/miplearn/solvers/gurobi.py @@ -13,13 +13,25 @@ logger = logging.getLogger(__name__) class GurobiSolver(InternalSolver): - def __init__(self, params=None): + def __init__(self, + params=None, + lazy_cb_frequency=1, + ): + """ + An InternalSolver backed by Gurobi's Python API (without Pyomo). + + Parameters + ---------- + params + Parameters to pass to Gurobi. For example, params={"MIPGap": 1e-3} + sets the gap tolerance to 1e-3. + lazy_cb_frequency + If 1, calls lazy constraint callbacks whenever an integer solution + is found. If 2, calls it also at every node, after solving the + LP relaxation of that node. + """ if params is None: params = {} - # params = { - # "LazyConstraints": 1, - # "PreCrush": 1, - # } from gurobipy import GRB self.GRB = GRB self.instance = None @@ -27,8 +39,16 @@ class GurobiSolver(InternalSolver): self.params = params self._all_vars = None self._bin_vars = None + self.cb_where = None + assert lazy_cb_frequency in [1, 2] + if lazy_cb_frequency == 1: + self.lazy_cb_where = [self.GRB.Callback.MIPSOL] + else: + self.lazy_cb_where = [self.GRB.Callback.MIPSOL, + self.GRB.Callback.MIPNODE] def set_instance(self, instance, model=None): + self._raise_if_callback() if model is None: model = instance.to_model() self.instance = instance @@ -36,6 +56,10 @@ class GurobiSolver(InternalSolver): self.model.update() self._update_vars() + def _raise_if_callback(self): + if self.cb_where is not None: + raise Exception("method cannot be called from a callback") + def _update_vars(self): self._all_vars = {} self._bin_vars = {} @@ -63,6 +87,7 @@ class GurobiSolver(InternalSolver): self.model.setParam(name, value) def solve_lp(self, tee=False): + self._raise_if_callback() self._apply_params() streams = [StringIO()] if tee: @@ -83,7 +108,24 @@ class GurobiSolver(InternalSolver): "Log": log } - def solve(self, tee=False, iteration_cb=None): + def solve(self, + tee=False, + iteration_cb=None, + lazy_cb=None): + self._raise_if_callback() + + def cb_wrapper(cb_model, cb_where): + try: + self.cb_where = cb_where + if cb_where in self.lazy_cb_where: + lazy_cb(self, self.model) + except: + logger.exception("callback error") + finally: + self.cb_where = None + + if lazy_cb: + self.params["LazyConstraints"] = 1 self._apply_params() total_wallclock_time = 0 total_nodes = 0 @@ -95,7 +137,10 @@ class GurobiSolver(InternalSolver): while True: logger.debug("Solving MIP...") with RedirectOutput(streams): - self.model.optimize() + if lazy_cb is None: + self.model.optimize() + else: + self.model.optimize(cb_wrapper) total_wallclock_time += self.model.runtime total_nodes += int(self.model.nodeCount) should_repeat = iteration_cb() @@ -114,6 +159,8 @@ class GurobiSolver(InternalSolver): } def get_solution(self): + self._raise_if_callback() + solution = {} for (varname, vardict) in self._all_vars.items(): solution[varname] = {} @@ -121,7 +168,22 @@ class GurobiSolver(InternalSolver): solution[varname][idx] = var.x return solution + def get_value(self, var_name, index): + var = self._all_vars[var_name][index] + return self._get_value(var) + + def _get_value(self, var): + if self.cb_where == self.GRB.Callback.MIPSOL: + return self.model.cbGetSolution(var) + elif self.cb_where == self.GRB.Callback.MIPNODE: + return self.model.cbGetNodeRel(var) + elif self.cb_where is None: + return var.x + else: + raise Exception("get_value cannot be called from cb_where=%s" % self.cb_where) + def get_variables(self): + self._raise_if_callback() variables = {} for (varname, vardict) in self._all_vars.items(): variables[varname] = [] @@ -132,12 +194,18 @@ class GurobiSolver(InternalSolver): def add_constraint(self, constraint, name=""): if type(constraint) is tuple: lhs, sense, rhs, name = constraint - logger.debug(lhs, sense, rhs) - self.model.addConstr(lhs, sense, rhs, name) + if self.cb_where in [self.GRB.Callback.MIPSOL, self.GRB.Callback.MIPNODE]: + self.model.cbLazy(lhs, sense, rhs) + else: + self.model.addConstr(lhs, sense, rhs, name) else: - self.model.addConstr(constraint, name=name) + if self.cb_where in [self.GRB.Callback.MIPSOL, self.GRB.Callback.MIPNODE]: + self.model.cbLazy(constraint) + else: + self.model.addConstr(constraint, name=name) def set_warm_start(self, solution): + self._raise_if_callback() count_fixed, count_total = 0, 0 for (varname, vardict) in solution.items(): for (idx, value) in vardict.items(): @@ -149,11 +217,13 @@ class GurobiSolver(InternalSolver): (count_fixed, count_total)) def clear_warm_start(self): + self._raise_if_callback() for (varname, vardict) in self._all_vars: for (idx, var) in vardict.items(): var[idx].start = self.GRB.UNDEFINED def fix(self, solution): + self._raise_if_callback() for (varname, vardict) in solution.items(): for (idx, value) in vardict.items(): if value is None: @@ -164,10 +234,12 @@ class GurobiSolver(InternalSolver): var.ub = value def get_constraint_ids(self): + self._raise_if_callback() self.model.update() return [c.ConstrName for c in self.model.getConstrs()] def extract_constraint(self, cid): + self._raise_if_callback() constr = self.model.getConstrByName(cid) cobj = (self.model.getRow(constr), constr.sense, @@ -178,29 +250,41 @@ class GurobiSolver(InternalSolver): def is_constraint_satisfied(self, cobj, tol=1e-5): lhs, sense, rhs, name = cobj - lhs_value = lhs.getValue() + if self.cb_where is not None: + lhs_value = lhs.getConstant() + for i in range(lhs.size()): + var = lhs.getVar(i) + coeff = lhs.getCoeff(i) + lhs_value += self._get_value(var) * coeff + else: + lhs_value = lhs.getValue() if sense == "<": return lhs_value <= rhs + tol elif sense == ">": return lhs_value >= rhs - tol elif sense == "=": - return abs(rhs - lhs_value) < tol + return abs(rhs - lhs_value) < abs(tol) else: raise Exception("Unknown sense: %s" % sense) def set_branching_priorities(self, priorities): + self._raise_if_callback() logger.warning("set_branching_priorities not implemented") def set_threads(self, threads): + self._raise_if_callback() self.params["Threads"] = threads def set_time_limit(self, time_limit): + self._raise_if_callback() self.params["TimeLimit"] = time_limit def set_node_limit(self, node_limit): + self._raise_if_callback() self.params["NodeLimit"] = node_limit def set_gap_tolerance(self, gap_tolerance): + self._raise_if_callback() self.params["MIPGap"] = gap_tolerance def _extract_warm_start_value(self, log): diff --git a/miplearn/solvers/internal.py b/miplearn/solvers/internal.py index 85bad83..9c9ee4a 100644 --- a/miplearn/solvers/internal.py +++ b/miplearn/solvers/internal.py @@ -119,7 +119,7 @@ class InternalSolver(ABC): pass @abstractmethod - def solve(self, tee=False, iteration_cb=None): + def solve(self, tee=False, iteration_cb=None, lazy_cb=None): """ Solves the currently loaded instance. After this method finishes, the best solution found can be retrieved by calling `get_solution`. @@ -132,7 +132,15 @@ class InternalSolver(ABC): instead, InternalSolver enters a loop, where `solve` and `iteration_cb` are called alternatively. To stop the loop, `iteration_cb` should return False. Any other result causes the solver to loop again. - tee: bool + lazy_cb: (internal_solver, model) -> None + This function is called whenever the solver finds a new candidate + solution and can be used to add lazy constraints to the model. Only + two operations within the callback are allowed: + - Querying the value of a variable, through `get_value(var, idx)` + - Querying if a constraint is satisfied, through `is_constraint_satisfied(cobj)` + - Adding a new constraint to the problem, through `add_constraint` + Additional operations may be allowed by specific subclasses. + tee: Bool If true, prints the solver log to the screen. Returns @@ -144,6 +152,13 @@ class InternalSolver(ABC): """ pass + @abstractmethod + def get_value(self, var_name, index): + """ + Returns the current value of a decision variable. + """ + pass + @abstractmethod def get_constraint_ids(self): """ diff --git a/miplearn/solvers/learning.py b/miplearn/solvers/learning.py index d8ae4b4..cf4eeef 100644 --- a/miplearn/solvers/learning.py +++ b/miplearn/solvers/learning.py @@ -31,34 +31,61 @@ def _parallel_solve(instance_idx): "Solution": instance.solution, "LP solution": instance.lp_solution, "Violated lazy constraints": instance.found_violated_lazy_constraints, - "Violated user cuts": instance.found_violated_user_cuts, + #"Violated user cuts": instance.found_violated_user_cuts, } class LearningSolver: - """ - Mixed-Integer Linear Programming (MIP) solver that extracts information - from previous runs, using Machine Learning methods, to accelerate the - solution of new (yet unseen) instances. - - Parameters - ---------- - solve_lp_first: bool - If true, solve LP relaxation first, then solve original MILP. This - option should be activated if the LP relaxation is not very - expensive to solve and if it provides good hints for the integer - solution. - """ - def __init__(self, components=None, - gap_tolerance=None, + gap_tolerance=1e-4, mode="exact", solver="gurobi", threads=None, time_limit=None, node_limit=None, - solve_lp_first=True): + solve_lp_first=True, + use_lazy_cb=False): + """ + Mixed-Integer Linear Programming (MIP) solver that extracts information + from previous runs and uses Machine Learning methods to accelerate the + solution of new (yet unseen) instances. + + Parameters + ---------- + components + Set of components in the solver. By default, includes: + - ObjectiveValueComponent + - PrimalSolutionComponent + - DynamicLazyConstraintsComponent + - UserCutsComponent + gap_tolerance + Relative MIP gap tolerance. By default, 1e-4. + mode + If "exact", solves problem to optimality, keeping all optimality + guarantees provided by the MIP solver. If "heuristic", uses machine + learning more agressively, and may return suboptimal solutions. + solver + The internal MIP solver to use. Can be either "cplex", "gurobi", a + solver class such as GurobiSolver, or a solver instance such as + GurobiSolver(). + threads + Maximum number of threads to use. If None, uses solver default. + time_limit + Maximum running time in seconds. If None, uses solver default. + node_limit + Maximum number of branch-and-bound nodes to explore. If None, uses + solver default. + use_lazy_cb + If True, uses lazy callbacks to enforce lazy constraints, instead of + a simple solver loop. This functionality may not supported by + all internal MIP solvers. + solve_lp_first: bool + If true, solve LP relaxation first, then solve original MILP. This + option should be activated if the LP relaxation is not very + expensive to solve and if it provides good hints for the integer + solution. + """ self.components = {} self.mode = mode self.internal_solver = None @@ -69,6 +96,7 @@ class LearningSolver: self.tee = False self.node_limit = node_limit self.solve_lp_first = solve_lp_first + self.use_lazy_cb = use_lazy_cb if components is not None: for comp in components: @@ -122,14 +150,12 @@ class LearningSolver: - instance.lower_bound - instance.upper_bound - instance.solution - - instance.found_violated_lazy_constraints - instance.solver_log - Additional solver components may set additional properties. Please see their documentation for more details. - If `solve_lp_first` is False, the properties lp_solution and lp_value - will be set to dummy values. + If `solver.solve_lp_first` is False, the properties lp_solution and + lp_value will be set to dummy values. Parameters ---------- @@ -175,13 +201,23 @@ class LearningSolver: def iteration_cb(): should_repeat = False - for component in self.components.values(): - if component.after_iteration(self, instance, model): + for comp in self.components.values(): + if comp.after_iteration(self, instance, model): should_repeat = True return should_repeat + def lazy_cb_wrapper(cb_solver, cb_model): + for comp in self.components.values(): + comp.on_lazy_callback(self, instance, model) + + lazy_cb = None + if self.use_lazy_cb: + lazy_cb = lazy_cb_wrapper + logger.info("Solving MILP...") - results = self.internal_solver.solve(tee=tee, iteration_cb=iteration_cb) + results = self.internal_solver.solve(tee=tee, + iteration_cb=iteration_cb, + lazy_cb=lazy_cb) results["LP value"] = instance.lp_value # Read MIP solution and bounds @@ -217,7 +253,7 @@ class LearningSolver: instances[idx].lower_bound = r["Results"]["Lower bound"] instances[idx].upper_bound = r["Results"]["Upper bound"] instances[idx].found_violated_lazy_constraints = r["Violated lazy constraints"] - instances[idx].found_violated_user_cuts = r["Violated user cuts"] + #instances[idx].found_violated_user_cuts = r["Violated user cuts"] instances[idx].solver_log = r["Results"]["Log"] return results diff --git a/miplearn/solvers/pyomo/base.py b/miplearn/solvers/pyomo/base.py index 89bb48e..b509985 100644 --- a/miplearn/solvers/pyomo/base.py +++ b/miplearn/solvers/pyomo/base.py @@ -57,6 +57,10 @@ class BasePyomoSolver(InternalSolver): solution[str(var)][index] = var[index].value return solution + def get_value(self, var_name, index): + var = self._varname_to_var[var_name] + return var[index].value + def get_variables(self): variables = {} for var in self.model.component_objects(Var): @@ -137,7 +141,12 @@ class BasePyomoSolver(InternalSolver): self._pyomo_solver.add_constraint(constraint) self._update_constrs() - def solve(self, tee=False, iteration_cb=None): + def solve(self, + tee=False, + iteration_cb=None, + lazy_cb=None): + if lazy_cb is not None: + raise Exception("lazy callback not supported") total_wallclock_time = 0 streams = [StringIO()] if tee: diff --git a/miplearn/solvers/tests/__init__.py b/miplearn/solvers/tests/__init__.py index 87a453a..870fa8f 100644 --- a/miplearn/solvers/tests/__init__.py +++ b/miplearn/solvers/tests/__init__.py @@ -7,13 +7,13 @@ from miplearn.problems.knapsack import KnapsackInstance, GurobiKnapsackInstance def _get_instance(solver): - if issubclass(solver, BasePyomoSolver): + if issubclass(solver, BasePyomoSolver) or isinstance(solver, BasePyomoSolver): return KnapsackInstance( weights=[23., 26., 20., 18.], prices=[505., 352., 458., 220.], capacity=67., ) - if issubclass(solver, GurobiSolver): + if issubclass(solver, GurobiSolver) or isinstance(solver, GurobiSolver): return GurobiKnapsackInstance( weights=[23., 26., 20., 18.], prices=[505., 352., 458., 220.], diff --git a/miplearn/solvers/tests/test_lazy_cb.py b/miplearn/solvers/tests/test_lazy_cb.py new file mode 100644 index 0000000..03e9861 --- /dev/null +++ b/miplearn/solvers/tests/test_lazy_cb.py @@ -0,0 +1,27 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +import logging + +from . import _get_instance +from ... import GurobiSolver + +logger = logging.getLogger(__name__) + + +def test_lazy_cb(): + solver = GurobiSolver() + instance = _get_instance(solver) + model = instance.to_model() + + def lazy_cb(cb_solver, cb_model): + logger.info("x[0] = %.f" % cb_solver.get_value("x", 0)) + cobj = (cb_model.getVarByName("x[0]") * 1.0, "<", 0.0, "cut") + if not cb_solver.is_constraint_satisfied(cobj): + cb_solver.add_constraint(cobj) + + solver.set_instance(instance, model) + solver.solve(lazy_cb=lazy_cb) + solution = solver.get_solution() + assert solution["x"][0] == 0.0