@@ -190,8 +186,43 @@

alt

-

Multidimensional 0-1 Knapsack Problem

+

Traveling Salesman Problem

Problem definition

+

Given a list of cities and the distance between each pair of cities, the problem asks for the +shortest route starting at the first city, visiting each other city exactly once, then returning +to the first city. This problem is a generalization of the Hamiltonian path problem, one of Karp's +21 NP-complete problems.

+

Random problem generator

+

The class TravelingSalesmanGenerator can be used to generate random instances of this +problem. Initially, the generator creates $n$ cities $(x_1,y_1),\ldots,(x_n,y_n) \in \mathbb{R}^2$, +where $n, x_i$ and $y_i$ are sampled independently from the provided probability distributions n, +x and y. For each pair of cities $(i,j)$, the distance $d_{i,j}$ between them is set to: + +where $\gamma_{i,j}$ is sampled from the distribution gamma.

+

If fix_cities=True is provided, the list of cities is kept the same for all generated instances. +The $gamma$ values, and therefore also the distances, are still different.

+

By default, all distances $d_{i,j}$ are rounded to the nearest integer. If round=False +is provided, this rounding will be disabled.

+

Challenge A

+ +
TravelingSalesmanGenerator(x=uniform(loc=0.0, scale=1000.0),
+                           y=uniform(loc=0.0, scale=1000.0),
+                           n=randint(low=350, high=351),
+                           gamma=uniform(loc=0.95, scale=0.1),
+                           fix_cities=True,
+                           round=True,
+                          )
+
+ +

alt

+

Multidimensional 0-1 Knapsack Problem

+

Problem definition

Given a set of $n$ items and $m$ types of resources (also called knapsacks), the problem is to find a subset of items that maximizes profit without consuming more resources than it is available. More precisely, the problem is:

- - - @@ -182,12 +173,12 @@ for instance in all_instances:

The first method is used by LearningSolver to construct a concrete Pyomo model, which will be provided to the internal MIP solver. The user should keep a reference to this Pyomo model, in order to retrieve, for example, the optimal variable values.

The second and third methods provide an encoding of the instance, which can be used by the ML models to make predictions. In the knapsack problem, for example, an implementation may decide to provide as instance features the average weights, average prices, number of items and the size of the knapsack. The weight and the price of each individual item could be provided as variable features. See miplearn/problems/knapsack.py for a concrete example.

-

An optional method which can be implemented is instance.get_variable_category(var, index), which returns a category (a string, an integer or any hashable type) for each decision variable. If two variables have the same category, LearningSolver will use the same internal ML model to predict the values of both variables. By default, all variables belong to the "default" category, and therefore only one ML model is used for all variables. If the returned category is None, ML predictors will ignore the variable.

-

It is not necessary to have a one-to-one correspondence between features and problem instances. One important (and deliberate) limitation of MIPLearn, however, is that get_instance_features() must always return arrays of same length for all relevant instances of the problem. Similarly, get_variable_features(var, index) must also always return arrays of same length for all variables in each category. It is up to the user to decide how to encode variable-length characteristics of the problem into fixed-length vectors. In graph problems, for example, graph embeddings can be used to reduce the (variable-length) lists of nodes and edges into a fixed-length structure that still preserves some properties of the graph. Different instance encodings may have significant impact on performance.

+

An optional method which can be implemented is instance.get_variable_category(var_name, index), which returns a category (a string, an integer or any hashable type) for each decision variable. If two variables have the same category, LearningSolver will use the same internal ML model to predict the values of both variables. By default, all variables belong to the "default" category, and therefore only one ML model is used for all variables. If the returned category is None, ML predictors will ignore the variable.

+

It is not necessary to have a one-to-one correspondence between features and problem instances. One important (and deliberate) limitation of MIPLearn, however, is that get_instance_features() must always return arrays of same length for all relevant instances of the problem. Similarly, get_variable_features(var_name, index) must also always return arrays of same length for all variables in each category. It is up to the user to decide how to encode variable-length characteristics of the problem into fixed-length vectors. In graph problems, for example, graph embeddings can be used to reduce the (variable-length) lists of nodes and edges into a fixed-length structure that still preserves some properties of the graph. Different instance encodings may have significant impact on performance.

Obtaining heuristic solutions

By default, LearningSolver uses Machine Learning to accelerate the MIP solution process, while maintaining all optimality guarantees provided by the MIP solver. In the default mode of operation, for example, predicted optimal solutions are used only as MIP starts.

For more significant performance benefits, LearningSolver can also be configured to place additional trust in the Machine Learning predictors, by using the mode="heuristic" constructor argument. When operating in this mode, if a ML model is statistically shown (through stratified k-fold cross validation) to have exceptionally high accuracy, the solver may decide to restrict the search space based on its predictions. The parts of the solution which the ML models cannot predict accurately will still be explored using traditional (branch-and-bound) methods. For particular applications, this mode has been shown to quickly produce optimal or near-optimal solutions (see references and benchmark results).

@@ -245,9 +236,10 @@ solver.solve(test_instance)

- Copyright © 2020, UChicago Argonne, LLC. All Rights Reserved.
 + Copyright © 2020, UChicago Argonne, LLC. All Rights Reserved.
- Documentation built with MkDocs.

+ Documentation built with MkDocs. +

diff --git a/miplearn/__init__.py b/miplearn/__init__.py index 8ee19d4..1c438aa 100644 --- a/miplearn/__init__.py +++ b/miplearn/__init__.py @@ -3,13 +3,13 @@ # Released under the modified BSD license. See COPYING.md for more details. from .extractors import (SolutionExtractor, - CombinedExtractor, InstanceFeaturesExtractor, ObjectiveValueExtractor, VariableFeaturesExtractor, ) from .components.component import Component from .components.objective import ObjectiveValueComponent +from .components.lazy import LazyConstraintsComponent from .components.primal import (PrimalSolutionComponent, AdaptivePredictor, ) diff --git a/miplearn/components/component.py b/miplearn/components/component.py index fe4ca6a..fba3bf1 100644 --- a/miplearn/components/component.py +++ b/miplearn/components/component.py @@ -18,10 +18,6 @@ class Component(ABC): def after_solve(self, solver, instance, model): pass - @abstractmethod - def merge(self, other): - pass - @abstractmethod def fit(self, training_instances): pass diff --git a/miplearn/components/lazy.py b/miplearn/components/lazy.py new file mode 100644 index 0000000..8200085 --- /dev/null +++ b/miplearn/components/lazy.py @@ -0,0 +1,57 @@ +# 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. + +from .component import Component +from ..extractors import * + +from abc import ABC, abstractmethod +from copy import deepcopy +import numpy as np +from sklearn.pipeline import make_pipeline +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler +from sklearn.model_selection import cross_val_score +from sklearn.metrics import roc_curve +from sklearn.neighbors import KNeighborsClassifier +from tqdm.auto import tqdm +import pyomo.environ as pe +import logging +logger = logging.getLogger(__name__) + + +class LazyConstraintsComponent(Component): + """ + A component that predicts which lazy constraints to enforce. + """ + + def __init__(self): + self.violations = set() + self.count = {} + self.n_samples = 0 + + def before_solve(self, solver, instance, model): + logger.info("Enforcing %d lazy constraints" % len(self.violations)) + for v in self.violations: + if self.count[v] < self.n_samples * 0.05: + continue + cut = instance.build_lazy_constraint(model, v) + solver.internal_solver.add_constraint(cut) + + def after_solve(self, solver, instance, model): + pass + + def fit(self, training_instances): + logger.debug("Fitting...") + self.n_samples = len(training_instances) + for instance in training_instances: + if not hasattr(instance, "found_violations"): + continue + for v in instance.found_violations: + self.violations.add(v) + if v not in self.count.keys(): + self.count[v] = 0 + self.count[v] += 1 + + def predict(self, instance, model=None): + return self.violations diff --git a/miplearn/components/objective.py b/miplearn/components/objective.py index ecf2f72..cf6d20e 100644 --- a/miplearn/components/objective.py +++ b/miplearn/components/objective.py @@ -30,16 +30,16 @@ class ObjectiveValueComponent(Component): def after_solve(self, solver, instance, model): pass - def merge(self, other): - pass - def fit(self, training_instances): + logger.debug("Extracting features...") features = InstanceFeaturesExtractor().extract(training_instances) ub = ObjectiveValueExtractor(kind="upper bound").extract(training_instances) lb = ObjectiveValueExtractor(kind="lower bound").extract(training_instances) self.ub_regressor = deepcopy(self.regressor_prototype) self.lb_regressor = deepcopy(self.regressor_prototype) + logger.debug("Fitting ub_regressor...") self.ub_regressor.fit(features, ub) + logger.debug("Fitting ub_regressor...") self.lb_regressor.fit(features, lb) def predict(self, instances): diff --git a/miplearn/components/primal.py b/miplearn/components/primal.py index 92d7542..3b21075 100644 --- a/miplearn/components/primal.py +++ b/miplearn/components/primal.py @@ -129,7 +129,7 @@ class PrimalSolutionComponent(Component): self.dynamic_thresholds = dynamic_thresholds def before_solve(self, solver, instance, model): - solution = self.predict(instance, model) + solution = self.predict(instance) if self.mode == "heuristic": solver.internal_solver.fix(solution) else: @@ -139,6 +139,7 @@ class PrimalSolutionComponent(Component): pass def fit(self, training_instances): + logger.debug("Extracting features...") features = VariableFeaturesExtractor().extract(training_instances) solutions = SolutionExtractor().extract(training_instances) @@ -180,12 +181,10 @@ class PrimalSolutionComponent(Component): self.thresholds[category, label] = thresholds[k] - def predict(self, instance, model=None): - if model is None: - model = instance.to_model() - x_test = VariableFeaturesExtractor().extract([instance], [model]) + def predict(self, instance): + x_test = VariableFeaturesExtractor().extract([instance]) solution = {} - var_split = Extractor.split_variables(instance, model) + var_split = Extractor.split_variables(instance) for category in var_split.keys(): for (i, (var, index)) in enumerate(var_split[category]): if var not in solution.keys(): @@ -200,6 +199,3 @@ class PrimalSolutionComponent(Component): if ws[i, 1] >= self.thresholds[category, label]: solution[var][index] = label return solution - - def merge(self, other_components): - pass diff --git a/miplearn/components/tests/test_primal.py b/miplearn/components/tests/test_primal.py index b35d5e7..3f9ec03 100644 --- a/miplearn/components/tests/test_primal.py +++ b/miplearn/components/tests/test_primal.py @@ -27,29 +27,7 @@ def test_predict(): instances, models = _get_instances() comp = PrimalSolutionComponent() comp.fit(instances) - solution = comp.predict(instances[0], models[0]) - assert models[0].x in solution.keys() + solution = comp.predict(instances[0]) + assert "x" in solution for idx in range(4): - assert idx in solution[models[0].x].keys() - -# def test_warm_start_save_load(): -# state_file = tempfile.NamedTemporaryFile(mode="r") -# solver = LearningSolver(components={"warm-start": WarmStartComponent()}) -# solver.parallel_solve(_get_instances(), n_jobs=2) -# solver.fit() -# comp = solver.components["warm-start"] -# assert comp.x_train["default"].shape == (8, 6) -# assert comp.y_train["default"].shape == (8, 2) -# assert ("default", 0) in comp.predictors.keys() -# assert ("default", 1) in comp.predictors.keys() -# solver.save_state(state_file.name) - -# solver.solve(_get_instances()[0]) - -# solver = LearningSolver(components={"warm-start": WarmStartComponent()}) -# solver.load_state(state_file.name) -# comp = solver.components["warm-start"] -# assert comp.x_train["default"].shape == (8, 6) -# assert comp.y_train["default"].shape == (8, 2) -# assert ("default", 0) in comp.predictors.keys() -# assert ("default", 1) in comp.predictors.keys() + assert idx in solution["x"] diff --git a/miplearn/extractors.py b/miplearn/extractors.py index f966250..d4d1d1d 100644 --- a/miplearn/extractors.py +++ b/miplearn/extractors.py @@ -5,6 +5,10 @@ import numpy as np from abc import ABC, abstractmethod from pyomo.core import Var +from tqdm.auto import tqdm, trange +from p_tqdm import p_map +import logging +logger = logging.getLogger(__name__) class Extractor(ABC): @@ -13,59 +17,39 @@ class Extractor(ABC): pass @staticmethod - def split_variables(instance, model): + def split_variables(instance): + assert hasattr(instance, "lp_solution") result = {} - for var in model.component_objects(Var): - for index in var: - category = instance.get_variable_category(var, index) + for var_name in instance.lp_solution: + for index in instance.lp_solution[var_name]: + category = instance.get_variable_category(var_name, index) if category is None: continue - if category not in result.keys(): + if category not in result: result[category] = [] - result[category] += [(var, index)] + result[category] += [(var_name, index)] return result - - @staticmethod - def merge(partial_results, vertical=False): - results = {} - all_categories = set() - for pr in partial_results: - all_categories |= pr.keys() - for category in all_categories: - results[category] = [] - for pr in partial_results: - if category in pr.keys(): - results[category] += [pr[category]] - if vertical: - results[category] = np.vstack(results[category]) - else: - results[category] = np.hstack(results[category]) - return results class VariableFeaturesExtractor(Extractor): - def extract(self, - instances, - models=None, - ): + def extract(self, instances): result = {} - if models is None: - models = [instance.to_model() for instance in instances] - for (index, instance) in enumerate(instances): - model = models[index] + for instance in tqdm(instances, + desc="Extract var features", + disable=len(instances) < 5): instance_features = instance.get_instance_features() - var_split = self.split_variables(instance, model) + var_split = self.split_variables(instance) for (category, var_index_pairs) in var_split.items(): - if category not in result.keys(): + if category not in result: result[category] = [] - for (var, index) in var_index_pairs: - result[category] += [np.hstack([ - instance_features, - instance.get_variable_features(var, index), - instance.lp_solution[str(var)][index], - ])] - for category in result.keys(): - result[category] = np.vstack(result[category]) + for (var_name, index) in var_index_pairs: + result[category] += [ + instance_features.tolist() + \ + instance.get_variable_features(var_name, index).tolist() + \ + [instance.lp_solution[var_name][index]] + ] + for category in result: + result[category] = np.array(result[category]) return result @@ -73,39 +57,29 @@ class SolutionExtractor(Extractor): def __init__(self, relaxation=False): self.relaxation = relaxation - def extract(self, instances, models=None): + def extract(self, instances): result = {} - if models is None: - models = [instance.to_model() for instance in instances] - for (index, instance) in enumerate(instances): - model = models[index] - var_split = self.split_variables(instance, model) + for instance in tqdm(instances, + desc="Extract solution", + disable=len(instances) < 5): + var_split = self.split_variables(instance) for (category, var_index_pairs) in var_split.items(): - if category not in result.keys(): + if category not in result: result[category] = [] - for (var, index) in var_index_pairs: + for (var_name, index) in var_index_pairs: if self.relaxation: - v = instance.lp_solution[str(var)][index] + v = instance.lp_solution[var_name][index] else: - v = instance.solution[str(var)][index] + v = instance.solution[var_name][index] if v is None: result[category] += [[0, 0]] else: result[category] += [[1 - v, v]] - for category in result.keys(): - result[category] = np.vstack(result[category]) + for category in result: + result[category] = np.array(result[category]) return result -class CombinedExtractor(Extractor): - def __init__(self, extractors): - self.extractors = extractors - - def extract(self, instances, models): - return self.merge([ex.extract(instances, models) - for ex in self.extractors]) - - class InstanceFeaturesExtractor(Extractor): def extract(self, instances, models=None): return np.vstack([ diff --git a/miplearn/instance.py b/miplearn/instance.py index a3fda45..40884fd 100644 --- a/miplearn/instance.py +++ b/miplearn/instance.py @@ -65,12 +65,50 @@ class Instance(ABC): def get_variable_category(self, var, index): """ - Returns a category (a string, an integer or any hashable type) for each decision variable. + Returns the category (a string, an integer or any hashable type) for each decision + variable. - If two variables have the same category, LearningSolver will use the same internal ML model - to predict the values of both variables. By default, all variables belong to the "default" - category, and therefore only one ML model is used for all variables. + If two variables have the same category, LearningSolver will use the same internal ML + model to predict the values of both variables. By default, all variables belong to the + "default" category, and therefore only one ML model is used for all variables. If the returned category is None, ML models will ignore the variable. """ return "default" + + def find_violations(self, model): + """ + Returns lazy constraint violations found for the current solution. + + After solving a model, LearningSolver will ask the instance to identify which lazy + constraints are violated by the current solution. For each identified violation, + LearningSolver will then call the build_lazy_constraint, add the generated Pyomo + constraint to the model, then resolve the problem. The process repeats until no further + lazy constraint violations are found. + + Each "violation" is simply a string, a tuple or any other hashable type which allows the + instance to identify unambiguously which lazy constraint should be generated. In the + Traveling Salesman Problem, for example, a subtour violation could be a frozen set + containing the cities in the subtour. + + For a concrete example, see TravelingSalesmanInstance. + """ + return [] + + def build_lazy_constraint(self, model, violation): + """ + Returns a Pyomo constraint which fixes a given violation. + + This method is typically called immediately after find_violations. The violation object + provided to this method is exactly the same object returned earlier by find_violations. + After some training, LearningSolver may decide to proactively build some lazy constraints + at the beginning of the optimization process, before a solution is even available. In this + case, build_lazy_constraints will be called without a corresponding call to + find_violations. + + The implementation should not directly add the constraint to the model. The constraint + will be added by LearningSolver after the method returns. + + For a concrete example, see TravelingSalesmanInstance. + """ + pass \ No newline at end of file diff --git a/miplearn/problems/tests/test_tsp.py b/miplearn/problems/tests/test_tsp.py new file mode 100644 index 0000000..d6182b5 --- /dev/null +++ b/miplearn/problems/tests/test_tsp.py @@ -0,0 +1,68 @@ +# 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. + +from miplearn import LearningSolver +from miplearn.problems.tsp import TravelingSalesmanGenerator, TravelingSalesmanInstance +import numpy as np +from numpy.linalg import norm +from scipy.spatial.distance import pdist, squareform +from scipy.stats import uniform, randint + + +def test_generator(): + instances = TravelingSalesmanGenerator(x=uniform(loc=0.0, scale=1000.0), + y=uniform(loc=0.0, scale=1000.0), + n=randint(low=100, high=101), + gamma=uniform(loc=0.95, scale=0.1), + fix_cities=True).generate(100) + assert len(instances) == 100 + assert instances[0].n_cities == 100 + assert norm(instances[0].distances - instances[0].distances.T) < 1e-6 + d = [instance.distances[0,1] for instance in instances] + assert np.std(d) > 0 + + +def test_instance(): + n_cities = 4 + distances = np.array([ + [0., 1., 2., 1.], + [1., 0., 1., 2.], + [2., 1., 0., 1.], + [1., 2., 1., 0.], + ]) + instance = TravelingSalesmanInstance(n_cities, distances) + solver = LearningSolver() + solver.solve(instance) + x = instance.solution["x"] + assert x[0,1] == 1.0 + assert x[0,2] == 0.0 + assert x[0,3] == 1.0 + assert x[1,2] == 1.0 + assert x[1,3] == 0.0 + assert x[2,3] == 1.0 + assert instance.lower_bound == 4.0 + assert instance.upper_bound == 4.0 + + +def test_subtour(): + n_cities = 6 + cities = np.array([ + [0., 0.], + [1., 0.], + [2., 0.], + [3., 0.], + [0., 1.], + [3., 1.], + ]) + distances = squareform(pdist(cities)) + instance = TravelingSalesmanInstance(n_cities, distances) + solver = LearningSolver() + solver.solve(instance) + x = instance.solution["x"] + assert x[0,1] == 1.0 + assert x[0,4] == 1.0 + assert x[1,2] == 1.0 + assert x[2,3] == 1.0 + assert x[3,5] == 1.0 + assert x[4,5] == 1.0 \ No newline at end of file diff --git a/miplearn/problems/tsp.py b/miplearn/problems/tsp.py new file mode 100644 index 0000000..e7fae64 --- /dev/null +++ b/miplearn/problems/tsp.py @@ -0,0 +1,169 @@ +# MIPLearn, an extensible framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2019-2020 Argonne National Laboratory. All rights reserved. +# Written by Alinson S. Xavier + +import numpy as np +import pyomo.environ as pe +from miplearn import Instance +from scipy.stats import uniform, randint +from scipy.spatial.distance import pdist, squareform +from scipy.stats.distributions import rv_frozen +import networkx as nx +import random + + +class ChallengeA: + def __init__(self, + seed=42, + n_training_instances=500, + n_test_instances=50, + ): + + np.random.seed(seed) + self.generator = TravelingSalesmanGenerator(x=uniform(loc=0.0, scale=1000.0), + y=uniform(loc=0.0, scale=1000.0), + n=randint(low=350, high=351), + gamma=uniform(loc=0.95, scale=0.1), + fix_cities=True, + round=True, + ) + + np.random.seed(seed + 1) + self.training_instances = self.generator.generate(n_training_instances) + + np.random.seed(seed + 2) + self.test_instances = self.generator.generate(n_test_instances) + + +class TravelingSalesmanGenerator: + """Random generator for the Traveling Salesman Problem.""" + + def __init__(self, + x=uniform(loc=0.0, scale=1000.0), + y=uniform(loc=0.0, scale=1000.0), + n=randint(low=100, high=101), + gamma=uniform(loc=1.0, scale=0.0), + fix_cities=True, + round=True, + ): + """Initializes the problem generator. + + Initially, the generator creates n cities (x_1,y_1),...,(x_n,y_n) where n, x_i and y_i are + sampled independently from the provided probability distributions `n`, `x` and `y`. For each + (unordered) pair of cities (i,j), the distance d[i,j] between them is set to: + + d[i,j] = gamma[i,j] \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2} + + where gamma is sampled from the provided probability distribution `gamma`. + + If fix_cities=True, the list of cities is kept the same for all generated instances. The + gamma values, and therefore also the distances, are still different. + + By default, all distances d[i,j] are rounded to the nearest integer. If `round=False` + is provided, this rounding will be disabled. + + Arguments + --------- + x: rv_continuous + Probability distribution for the x-coordinate of each city. + y: rv_continuous + Probability distribution for the y-coordinate of each city. + n: rv_discrete + Probability distribution for the number of cities. + fix_cities: bool + If False, cities will be resampled for every generated instance. Otherwise, list of + cities will be computed once, during the constructor. + round: bool + If True, distances are rounded to the nearest integer. + """ + assert isinstance(x, rv_frozen), "x should be a SciPy probability distribution" + assert isinstance(y, rv_frozen), "y should be a SciPy probability distribution" + assert isinstance(n, rv_frozen), "n should be a SciPy probability distribution" + assert isinstance(gamma, rv_frozen), "gamma should be a SciPy probability distribution" + self.x = x + self.y = y + self.n = n + self.gamma = gamma + self.round = round + + if fix_cities: + self.fixed_n, self.fixed_cities = self._generate_cities() + else: + self.fixed_n = None + self.fixed_cities = None + + def generate(self, n_samples): + def _sample(): + if self.fixed_cities is not None: + n, cities = self.fixed_n, self.fixed_cities + else: + n, cities = self._generate_cities() + distances = squareform(pdist(cities)) * self.gamma.rvs(size=(n, n)) + distances = np.tril(distances) + np.triu(distances.T, 1) + if self.round: + distances = distances.round() + return TravelingSalesmanInstance(n, distances) + return [_sample() for _ in range(n_samples)] + + def _generate_cities(self): + n = self.n.rvs() + cities = np.array([(self.x.rvs(), self.y.rvs()) for _ in range(n)]) + return n, cities + + +class TravelingSalesmanInstance(Instance): + """An instance ot the Traveling Salesman Problem. + + Given a list of cities and the distance between each pair of cities, the problem asks for the + shortest route starting at the first city, visiting each other city exactly once, then + returning to the first city. This problem is a generalization of the Hamiltonian path problem, + one of Karp's 21 NP-complete problems. + """ + + def __init__(self, n_cities, distances): + assert isinstance(distances, np.ndarray) + assert distances.shape == (n_cities, n_cities) + self.n_cities = n_cities + self.distances = distances + + def to_model(self): + model = pe.ConcreteModel() + model.edges = edges = [(i,j) + for i in range(self.n_cities) + for j in range(i+1, self.n_cities)] + model.x = pe.Var(edges, domain=pe.Binary) + model.obj = pe.Objective(expr=sum(model.x[i,j] * self.distances[i,j] + for (i,j) in edges), + sense=pe.minimize) + model.eq_degree = pe.ConstraintList() + model.eq_subtour = pe.ConstraintList() + for i in range(self.n_cities): + model.eq_degree.add(sum(model.x[min(i,j), max(i,j)] + for j in range(self.n_cities) if i != j) == 2) + return model + + def get_instance_features(self): + return np.array([1]) + + def get_variable_features(self, var_name, index): + return np.array([1]) + + def get_variable_category(self, var_name, index): + return index + + def find_violations(self, model): + selected_edges = [e for e in model.edges if model.x[e].value > 0.5] + graph = nx.Graph() + graph.add_edges_from(selected_edges) + components = [frozenset(c) for c in list(nx.connected_components(graph))] + violations = [] + for c in components: + if len(c) < self.n_cities: + violations += [c] + return violations + + def build_lazy_constraint(self, model, component): + cut_edges = [e for e in model.edges + if (e[0] in component and e[1] not in component) or + (e[0] not in component and e[1] in component)] + return model.eq_subtour.add(sum(model.x[e] for e in cut_edges) >= 2) diff --git a/miplearn/solvers.py b/miplearn/solvers.py index 7ff976c..9077639 100644 --- a/miplearn/solvers.py +++ b/miplearn/solvers.py @@ -2,45 +2,66 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -from . import ObjectiveValueComponent, PrimalSolutionComponent +from . import ObjectiveValueComponent, PrimalSolutionComponent, LazyConstraintsComponent import pyomo.environ as pe from pyomo.core import Var from copy import deepcopy import pickle from scipy.stats import randint from p_tqdm import p_map +import numpy as np import logging logger = logging.getLogger(__name__) +# Global memory for multiprocessing +SOLVER = [None] +INSTANCES = [None] + + +def _parallel_solve(instance_idx): + solver = deepcopy(SOLVER[0]) + instance = INSTANCES[0][instance_idx] + results = solver.solve(instance) + return { + "Results": results, + "Solution": instance.solution, + "LP solution": instance.lp_solution, + "LP value": instance.lp_value, + "Upper bound": instance.upper_bound, + "Lower bound": instance.lower_bound, + "Violations": instance.found_violations, + } + + class InternalSolver: def __init__(self): self.is_warm_start_available = False self.model = None - pass + self.var_name_to_var = {} def solve_lp(self, tee=False): + self.solver.set_instance(self.model) + # Relax domain from pyomo.core.base.set_types import Reals - original_domain = {} - for var in self.model.component_data_objects(Var): - original_domain[str(var)] = var.domain + original_domains = [] + for (idx, var) in enumerate(self.model.component_data_objects(Var)): + original_domains += [var.domain] lb, ub = var.bounds var.setlb(lb) var.setub(ub) var.domain = Reals + self.solver.update_var(var) # Solve LP relaxation - self.solver.set_instance(self.model) results = self.solver.solve(tee=tee) # Restore domains - for var in self.model.component_data_objects(Var): - var.domain = original_domain[str(var)] + for (idx, var) in enumerate(self.model.component_data_objects(Var)): + var.domain = original_domains[idx] + self.solver.update_var(var) - # Reload original model - self.solver.set_instance(self.model) - return { "Optimal value": results["Problem"][0]["Lower bound"], } @@ -58,36 +79,43 @@ class InternalSolver: solution[str(var)][index] = var[index].value return solution - def set_warm_start(self, ws): + def set_warm_start(self, solution): self.is_warm_start_available = True self.clear_values() count_total, count_fixed = 0, 0 - for var in ws.keys(): - for index in var: + for var_name in solution: + var = self.var_name_to_var[var_name] + for index in solution[var_name]: count_total += 1 - var[index].value = ws[var][index] - if ws[var][index] is not None: + var[index].value = solution[var_name][index] + if solution[var_name][index] is not None: count_fixed += 1 logger.info("Setting start values for %d variables (out of %d)" % (count_fixed, count_total)) - def set_model(self, model): self.model = model self.solver.set_instance(model) + self.var_name_to_var = {} + for var in model.component_objects(Var): + self.var_name_to_var[var.name] = var - def fix(self, ws): + def fix(self, solution): count_total, count_fixed = 0, 0 - for var in ws.keys(): - for index in var: + for var_name in solution: + for index in solution[var_name]: + var = self.var_name_to_var[var_name] count_total += 1 - if ws[var][index] is None: + if solution[var_name][index] is None: continue count_fixed += 1 - var[index].fix(ws[var][index]) + var[index].fix(solution[var_name][index]) self.solver.update_var(var[index]) logger.info("Fixing values for %d variables (out of %d)" % (count_fixed, count_total)) + + def add_constraint(self, cut): + self.solver.add_constraint(cut) class GurobiSolver(InternalSolver): @@ -198,6 +226,7 @@ class LearningSolver: self.components = { "ObjectiveValue": ObjectiveValueComponent(), "PrimalSolution": PrimalSolutionComponent(), + "LazyConstraints": LazyConstraintsComponent(), } assert self.mode in ["exact", "heuristic"] @@ -231,27 +260,44 @@ class LearningSolver: self.internal_solver = self._create_internal_solver() self.internal_solver.set_model(model) - # Solve LP relaxation + logger.debug("Solving LP relaxation...") results = self.internal_solver.solve_lp(tee=tee) instance.lp_solution = self.internal_solver.get_solution() instance.lp_value = results["Optimal value"] - # Invoke before_solve callbacks + logger.debug("Running before_solve callbacks...") for component in self.components.values(): component.before_solve(self, instance, model) if relaxation_only: return results - # Solver original MIP - results = self.internal_solver.solve(tee=tee) + total_wallclock_time = 0 + instance.found_violations = [] + while True: + logger.debug("Solving MIP...") + results = self.internal_solver.solve(tee=tee) + logger.debug(" %.2f s" % results["Wallclock time"]) + total_wallclock_time += results["Wallclock time"] + if not hasattr(instance, "find_violations"): + break + logger.debug("Finding violated constraints...") + violations = instance.find_violations(model) + if len(violations) == 0: + break + instance.found_violations += violations + logger.debug(" %d violations found" % len(violations)) + for v in violations: + cut = instance.build_lazy_constraint(model, v) + self.internal_solver.add_constraint(cut) + results["Wallclock time"] = total_wallclock_time # Read MIP solution and bounds instance.lower_bound = results["Lower bound"] instance.upper_bound = results["Upper bound"] instance.solution = self.internal_solver.get_solution() - # Invoke after_solve callbacks + logger.debug("Calling after_solve callbacks...") for component in self.components.values(): component.after_solve(self, instance, model) @@ -266,40 +312,23 @@ class LearningSolver: label="Solve", collect_training_data=True, ): - self.internal_solver = None - def _process(instance): - solver = deepcopy(self) - results = solver.solve(instance) - solver.internal_solver = None - if not collect_training_data: - solver.components = {} - return { - "Solver": solver, - "Results": results, - "Solution": instance.solution, - "LP solution": instance.lp_solution, - "LP value": instance.lp_value, - "Upper bound": instance.upper_bound, - "Lower bound": instance.lower_bound, - } + self.internal_solver = None + SOLVER[0] = self + INSTANCES[0] = instances + p_map_results = p_map(_parallel_solve, + list(range(len(instances))), + num_cpus=n_jobs, + desc=label) - p_map_results = p_map(_process, instances, num_cpus=n_jobs, desc=label) - subsolvers = [p["Solver"] for p in p_map_results] results = [p["Results"] for p in p_map_results] - for (idx, r) in enumerate(p_map_results): instances[idx].solution = r["Solution"] instances[idx].lp_solution = r["LP solution"] instances[idx].lp_value = r["LP value"] instances[idx].lower_bound = r["Lower bound"] instances[idx].upper_bound = r["Upper bound"] - - for (name, component) in self.components.items(): - subcomponents = [subsolver.components[name] - for subsolver in subsolvers - if name in subsolver.components.keys()] - self.components[name].merge(subcomponents) + instances[idx].found_violations = r["Violations"] return results @@ -310,21 +339,3 @@ class LearningSolver: return for component in self.components.values(): component.fit(training_instances) - - def save_state(self, filename): - with open(filename, "wb") as file: - pickle.dump({ - "version": 2, - "components": self.components, - }, file) - - def load_state(self, filename): - with open(filename, "rb") as file: - data = pickle.load(file) - assert data["version"] == 2 - for (component_name, component) in data["components"].items(): - if component_name not in self.components.keys(): - continue - else: - self.components[component_name].merge([component]) - diff --git a/miplearn/tests/test_benchmark.py b/miplearn/tests/test_benchmark.py index 094eee4..523ba61 100644 --- a/miplearn/tests/test_benchmark.py +++ b/miplearn/tests/test_benchmark.py @@ -18,8 +18,6 @@ def test_benchmark(): # Training phase... training_solver = LearningSolver() training_solver.parallel_solve(train_instances, n_jobs=10) - training_solver.fit() - training_solver.save_state("data.bin") # Test phase... test_solvers = { @@ -27,7 +25,7 @@ def test_benchmark(): "Strategy B": LearningSolver(), } benchmark = BenchmarkRunner(test_solvers) - benchmark.load_state("data.bin") + benchmark.fit(train_instances) benchmark.parallel_solve(test_instances, n_jobs=2, n_trials=2) assert benchmark.raw_results().values.shape == (12,13) diff --git a/miplearn/tests/test_extractors.py b/miplearn/tests/test_extractors.py index eb57df7..a51b386 100644 --- a/miplearn/tests/test_extractors.py +++ b/miplearn/tests/test_extractors.py @@ -5,7 +5,6 @@ from miplearn.problems.knapsack import KnapsackInstance from miplearn import (LearningSolver, SolutionExtractor, - CombinedExtractor, InstanceFeaturesExtractor, VariableFeaturesExtractor, ) @@ -33,7 +32,7 @@ def _get_instances(): def test_solution_extractor(): instances, models = _get_instances() - features = SolutionExtractor().extract(instances, models) + features = SolutionExtractor().extract(instances) assert isinstance(features, dict) assert "default" in features.keys() assert isinstance(features["default"], np.ndarray) @@ -48,17 +47,6 @@ def test_solution_extractor(): ] -def test_combined_extractor(): - instances, models = _get_instances() - extractor = CombinedExtractor(extractors=[VariableFeaturesExtractor(), - SolutionExtractor()]) - features = extractor.extract(instances, models) - assert isinstance(features, dict) - assert "default" in features.keys() - assert isinstance(features["default"], np.ndarray) - assert features["default"].shape == (6, 7) - - def test_instance_features_extractor(): instances, models = _get_instances() features = InstanceFeaturesExtractor().extract(instances) diff --git a/miplearn/tests/test_solver.py b/miplearn/tests/test_solver.py index c8b662e..24d879b 100644 --- a/miplearn/tests/test_solver.py +++ b/miplearn/tests/test_solver.py @@ -41,29 +41,6 @@ def test_solver(): solver.fit() solver.solve(instance) - -# def test_solve_save_load_state(): -# instance = _get_instance() -# components_before = { -# "warm-start": WarmStartComponent(), -# } -# solver = LearningSolver(components=components_before) -# solver.solve(instance) -# solver.fit() -# solver.save_state("/tmp/knapsack_train.bin") -# prev_x_train_len = len(solver.components["warm-start"].x_train) -# prev_y_train_len = len(solver.components["warm-start"].y_train) - -# components_after = { -# "warm-start": WarmStartComponent(), -# } -# solver = LearningSolver(components=components_after) -# solver.load_state("/tmp/knapsack_train.bin") -# assert len(solver.components.keys()) == 1 -# assert len(solver.components["warm-start"].x_train) == prev_x_train_len -# assert len(solver.components["warm-start"].y_train) == prev_y_train_len - - def test_parallel_solve(): instances = [_get_instance() for _ in range(10)] solver = LearningSolver()