From b5e602cdc16f9e9919f98857edea8463a49a3b1d Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Sat, 10 Apr 2021 17:24:03 -0500 Subject: [PATCH] get_constraints: Fetch slack and dual values --- miplearn/features.py | 11 +++-- miplearn/solvers/gurobi.py | 38 ++++++++++---- miplearn/solvers/internal.py | 8 --- miplearn/solvers/pyomo/base.py | 79 ++++++++++++++++++------------ miplearn/solvers/tests/__init__.py | 63 +++++++++++++++++++----- 5 files changed, 133 insertions(+), 66 deletions(-) diff --git a/miplearn/features.py b/miplearn/features.py index 6edffaf..1bbe957 100644 --- a/miplearn/features.py +++ b/miplearn/features.py @@ -43,12 +43,17 @@ class VariableFeatures: @dataclass class Constraint: - rhs: float = 0.0 + basis_status: Optional[str] = None + category: Optional[Hashable] = None + dual_value: Optional[float] = None + lazy: bool = False lhs: Dict[str, float] = lambda: {} # type: ignore + rhs: float = 0.0 + sa_rhs_down: Optional[float] = None + sa_rhs_up: Optional[float] = None sense: str = "<" + slack: Optional[float] = None user_features: Optional[List[float]] = None - lazy: bool = False - category: Hashable = None @dataclass diff --git a/miplearn/solvers/gurobi.py b/miplearn/solvers/gurobi.py index 0eaa149..da83837 100644 --- a/miplearn/solvers/gurobi.py +++ b/miplearn/solvers/gurobi.py @@ -5,6 +5,7 @@ import logging import re import sys from dataclasses import dataclass +from enum import Enum from io import StringIO from random import randint from typing import List, Any, Dict, Optional, Hashable @@ -68,6 +69,9 @@ class GurobiSolver(InternalSolver): self.lazy_cb_frequency = lazy_cb_frequency self._bin_vars: List["gurobipy.Var"] = [] self._varname_to_var: Dict[str, "gurobipy.Var"] = {} + self._dirty = True + self._has_lp_solution = False + self._has_mip_solution = False if self.lazy_cb_frequency == 1: self.lazy_cb_where = [self.gp.GRB.Callback.MIPSOL] @@ -136,9 +140,12 @@ class GurobiSolver(InternalSolver): var.ub = 1.0 with _RedirectOutput(streams): self.model.optimize() + self._dirty = False for var in self._bin_vars: var.vtype = self.gp.GRB.BINARY log = streams[0].getvalue() + self._has_lp_solution = self.model.solCount > 0 + self._has_mip_solution = False opt_value = None if not self.is_infeasible(): opt_value = self.model.objVal @@ -191,6 +198,7 @@ class GurobiSolver(InternalSolver): while True: with _RedirectOutput(streams): self.model.optimize(cb_wrapper) + self._dirty = False if len(callback_exceptions) > 0: raise callback_exceptions[0] total_wallclock_time += self.model.runtime @@ -198,6 +206,8 @@ class GurobiSolver(InternalSolver): should_repeat = iteration_cb() if not should_repeat: break + self._has_lp_solution = False + self._has_mip_solution = self.model.solCount > 0 # Fetch results and stats log = streams[0].getvalue() @@ -306,6 +316,9 @@ class GurobiSolver(InternalSolver): self.model.addConstr(lhs <= constr.rhs, name=name) else: self.model.addConstr(lhs >= constr.rhs, name=name) + self._dirty = True + self._has_lp_solution = False + self._has_mip_solution = False @overrides def remove_constraint(self, name: str) -> None: @@ -341,12 +354,6 @@ class GurobiSolver(InternalSolver): var.lb = value var.ub = value - @overrides - def get_inequality_slacks(self) -> Dict[str, float]: - assert self.model is not None - ineqs = [c for c in self.model.getConstrs() if c.sense != "="] - return {c.ConstrName: c.Slack for c in ineqs} - @overrides def relax(self) -> None: assert self.model is not None @@ -414,7 +421,9 @@ class GurobiSolver(InternalSolver): def get_constraints(self) -> Dict[str, Constraint]: assert self.model is not None self._raise_if_callback() - self.model.update() + if self._dirty: + self.model.update() + self._dirty = False constraints: Dict[str, Constraint] = {} for c in self.model.getConstrs(): constr = self._parse_gurobi_constraint(c) @@ -422,13 +431,22 @@ class GurobiSolver(InternalSolver): constraints[c.constrName] = constr return constraints - def _parse_gurobi_constraint(self, c: Any) -> Constraint: + def _parse_gurobi_constraint(self, gp_constr: Any) -> Constraint: assert self.model is not None - expr = self.model.getRow(c) + expr = self.model.getRow(gp_constr) lhs: Dict[str, float] = {} for i in range(expr.size()): lhs[expr.getVar(i).varName] = expr.getCoeff(i) - return Constraint(rhs=c.rhs, lhs=lhs, sense=c.sense) + constr = Constraint( + rhs=gp_constr.rhs, + lhs=lhs, + sense=gp_constr.sense, + ) + if self._has_lp_solution: + constr.dual_value = gp_constr.pi + if self._has_lp_solution or self._has_mip_solution: + constr.slack = gp_constr.slack + return constr @overrides def are_callbacks_supported(self) -> bool: diff --git a/miplearn/solvers/internal.py b/miplearn/solvers/internal.py index 243b5db..0149549 100644 --- a/miplearn/solvers/internal.py +++ b/miplearn/solvers/internal.py @@ -184,14 +184,6 @@ class InternalSolver(ABC, EnforceOverrides): """ pass - @abstractmethod - def get_inequality_slacks(self) -> Dict[str, float]: - """ - Returns a dictionary mapping constraint name to the constraint slack - in the current solution. - """ - pass - @abstractmethod def is_infeasible(self) -> bool: """ diff --git a/miplearn/solvers/pyomo/base.py b/miplearn/solvers/pyomo/base.py index 1934bd8..2785985 100644 --- a/miplearn/solvers/pyomo/base.py +++ b/miplearn/solvers/pyomo/base.py @@ -12,7 +12,7 @@ import numpy as np import pyomo from overrides import overrides from pyomo import environ as pe -from pyomo.core import Var +from pyomo.core import Var, Suffix from pyomo.core.base import _GeneralVarData from pyomo.core.base.constraint import ConstraintList from pyomo.core.expr.numeric_expr import SumExpression, MonomialTermExpression @@ -61,6 +61,8 @@ class BasePyomoSolver(InternalSolver): self._varname_to_var: Dict[str, pe.Var] = {} self._cname_to_constr: Dict[str, pe.Constraint] = {} self._termination_condition: str = "" + self._has_lp_solution = False + self._has_mip_solution = False for (key, value) in params.items(): self._pyomo_solver.options[key] = value @@ -76,10 +78,14 @@ class BasePyomoSolver(InternalSolver): streams += [sys.stdout] with _RedirectOutput(streams): results = self._pyomo_solver.solve(tee=True) + self._termination_condition = results["Solver"][0]["Termination condition"] self._restore_integrality() opt_value = None + self._has_lp_solution = False + self._has_mip_solution = False if not self.is_infeasible(): opt_value = results["Problem"][0]["Lower bound"] + self._has_lp_solution = True return { "LP value": opt_value, "LP log": streams[0].getvalue(), @@ -122,7 +128,10 @@ class BasePyomoSolver(InternalSolver): ws_value = self._extract_warm_start_value(log) self._termination_condition = results["Solver"][0]["Termination condition"] lb, ub = None, None + self._has_mip_solution = False + self._has_lp_solution = False if not self.is_infeasible(): + self._has_mip_solution = True lb = results["Problem"][0]["Lower bound"] ub = results["Problem"][0]["Upper bound"] stats: MIPSolveStats = { @@ -185,6 +194,7 @@ class BasePyomoSolver(InternalSolver): self.instance = instance self.model = model self.model.extra_constraints = ConstraintList() + self.model.dual = Suffix(direction=Suffix.IMPORT) self._pyomo_solver.set_instance(model) self._update_obj() self._update_vars() @@ -256,6 +266,9 @@ class BasePyomoSolver(InternalSolver): self._cname_to_constr[name] = cl else: self._pyomo_solver.add_constraint(constr) + self._termination_condition = "" + self._has_lp_solution = False + self._has_mip_solution = False @overrides def remove_constraint(self, name: str) -> None: @@ -321,22 +334,14 @@ class BasePyomoSolver(InternalSolver): var.domain = pyomo.core.base.set_types.Reals self._pyomo_solver.update_var(var) - @overrides - def get_inequality_slacks(self) -> Dict[str, float]: - result: Dict[str, float] = {} - for (cname, cobj) in self._cname_to_constr.items(): - if cobj.equality: - continue - result[cname] = cobj.slack() - return result - @overrides def is_infeasible(self) -> bool: return self._termination_condition == TerminationCondition.infeasible @overrides def get_dual(self, cid: str) -> float: - raise NotImplementedError() + constr = self._cname_to_constr[cid] + return self._pyomo_solver.dual[constr] @overrides def get_sense(self) -> str: @@ -376,45 +381,55 @@ class BasePyomoSolver(InternalSolver): return constraints - @staticmethod - def _parse_pyomo_constraint(c: pyomo.core.Constraint) -> Constraint: + def _parse_pyomo_constraint( + self, + pyomo_constr: pyomo.core.Constraint, + ) -> Constraint: + constr = Constraint() + # Extract RHS and sense - has_ub = c.has_ub() - has_lb = c.has_lb() + has_ub = pyomo_constr.has_ub() + has_lb = pyomo_constr.has_lb() assert ( - (not has_lb) or (not has_ub) or c.upper() == c.lower() + (not has_lb) or (not has_ub) or pyomo_constr.upper() == pyomo_constr.lower() ), "range constraints not supported" if has_lb: - sense = ">" - rhs = c.lower() + constr.sense = ">" + constr.rhs = pyomo_constr.lower() elif has_ub: - sense = "<" - rhs = c.upper() + constr.sense = "<" + constr.rhs = pyomo_constr.upper() else: - sense = "=" - rhs = c.upper() + constr.sense = "=" + constr.rhs = pyomo_constr.upper() # Extract LHS lhs = {} - if isinstance(c.body, SumExpression): - for term in c.body._args_: + if isinstance(pyomo_constr.body, SumExpression): + for term in pyomo_constr.body._args_: if isinstance(term, MonomialTermExpression): lhs[term._args_[1].name] = term._args_[0] elif isinstance(term, _GeneralVarData): lhs[term.name] = 1.0 else: raise Exception(f"Unknown term type: {term.__class__.__name__}") - elif isinstance(c.body, _GeneralVarData): - lhs[c.body.name] = 1.0 + elif isinstance(pyomo_constr.body, _GeneralVarData): + lhs[pyomo_constr.body.name] = 1.0 else: - raise Exception(f"Unknown expression type: {c.body.__class__.__name__}") + raise Exception( + f"Unknown expression type: {pyomo_constr.body.__class__.__name__}" + ) + constr.lhs = lhs + + # Extract solution attributes + if self._has_lp_solution: + constr.dual_value = self.model.dual[pyomo_constr] + + if self._has_mip_solution or self._has_lp_solution: + constr.slack = pyomo_constr.slack() # Build constraint - return Constraint( - lhs=lhs, - rhs=rhs, - sense=sense, - ) + return constr @overrides def are_callbacks_supported(self) -> bool: diff --git a/miplearn/solvers/tests/__init__.py b/miplearn/solvers/tests/__init__.py index 233c5c8..2637075 100644 --- a/miplearn/solvers/tests/__init__.py +++ b/miplearn/solvers/tests/__init__.py @@ -12,6 +12,14 @@ from miplearn.solvers.internal import InternalSolver # This file is in the main source folder, so that it can be called from Julia. +def _round_constraints(constraints): + for (cname, c) in constraints.items(): + for attr in ["slack", "dual_value"]: + if getattr(c, attr) is not None: + setattr(c, attr, round(getattr(c, attr), 6)) + return constraints + + def run_internal_solver_tests(solver: InternalSolver) -> None: run_basic_usage_tests(solver.clone()) run_warm_start_tests(solver.clone()) @@ -22,20 +30,38 @@ def run_internal_solver_tests(solver: InternalSolver) -> None: def run_basic_usage_tests(solver: InternalSolver) -> None: + # Create and set instance instance = solver.build_test_instance_knapsack() model = instance.to_model() solver.set_instance(instance, model) + + # Fetch variables (after-load) assert_equals( solver.get_variable_names(), ["x[0]", "x[1]", "x[2]", "x[3]"], ) + # Fetch constraints (after-load) + assert_equals( + _round_constraints(solver.get_constraints()), + { + "eq_capacity": Constraint( + lazy=False, + lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0}, + rhs=67.0, + sense="<", + ) + }, + ) + + # Solve linear programming relaxation lp_stats = solver.solve_lp() assert not solver.is_infeasible() assert lp_stats["LP value"] is not None assert_equals(round(lp_stats["LP value"], 3), 1287.923) assert len(lp_stats["LP log"]) > 100 + # Fetch variables (after-lp) solution = solver.get_solution() assert solution is not None assert solution["x[0]"] is not None @@ -47,6 +73,22 @@ def run_basic_usage_tests(solver: InternalSolver) -> None: assert_equals(round(solution["x[2]"], 3), 1.000) assert_equals(round(solution["x[3]"], 3), 0.000) + # Fetch constraints (after-lp) + assert_equals( + _round_constraints(solver.get_constraints()), + { + "eq_capacity": Constraint( + lazy=False, + lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0}, + rhs=67.0, + sense="<", + slack=0.0, + dual_value=13.538462, + ) + }, + ) + + # Solve MIP mip_stats = solver.solve( tee=True, iteration_cb=None, @@ -60,6 +102,7 @@ def run_basic_usage_tests(solver: InternalSolver) -> None: assert_equals(mip_stats["Sense"], "max") assert isinstance(mip_stats["Wallclock time"], float) + # Fetch variables (after-mip) solution = solver.get_solution() assert solution is not None assert solution["x[0]"] is not None @@ -71,13 +114,15 @@ def run_basic_usage_tests(solver: InternalSolver) -> None: assert_equals(solution["x[2]"], 1.0) assert_equals(solution["x[3]"], 1.0) + # Fetch constraints (after-mip) assert_equals( - solver.get_constraints(), + _round_constraints(solver.get_constraints()), { "eq_capacity": Constraint( lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0}, rhs=67.0, sense="<", + slack=6.0, ), }, ) @@ -86,10 +131,11 @@ def run_basic_usage_tests(solver: InternalSolver) -> None: cut = Constraint(lhs={"x[0]": 1.0}, sense="<", rhs=0.0) assert not solver.is_constraint_satisfied(cut) - # Add new constraint and verify that it is listed + # Add new constraint and verify that it is listed. Modifying the model should + # also clear the current solution. solver.add_constraint(cut, "cut") assert_equals( - solver.get_constraints(), + _round_constraints(solver.get_constraints()), { "eq_capacity": Constraint( lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0}, @@ -104,17 +150,11 @@ def run_basic_usage_tests(solver: InternalSolver) -> None: }, ) - # New constraint should affect the solution + # Re-solve MIP and verify that constraint affects the solution stats = solver.solve() assert_equals(stats["Lower bound"], 1030.0) assert solver.is_constraint_satisfied(cut) - # Verify slacks - assert_equals( - solver.get_inequality_slacks(), - {"cut": 0.0, "eq_capacity": 3.0}, - ) - # Remove the new constraint solver.remove_constraint("cut") @@ -122,9 +162,6 @@ def run_basic_usage_tests(solver: InternalSolver) -> None: stats = solver.solve() assert_equals(stats["Lower bound"], 1183.0) - # Constraint should not be satisfied by current solution - assert not solver.is_constraint_satisfied(cut) - def run_warm_start_tests(solver: InternalSolver) -> None: instance = solver.build_test_instance_knapsack()