diff --git a/pyomo/contrib/observer/__init__.py b/pyomo/contrib/observer/__init__.py new file mode 100644 index 00000000000..6eb9ea8b81d --- /dev/null +++ b/pyomo/contrib/observer/__init__.py @@ -0,0 +1,10 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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. +# ___________________________________________________________________________ diff --git a/pyomo/contrib/observer/component_collector.py b/pyomo/contrib/observer/component_collector.py new file mode 100644 index 00000000000..d52ec46086c --- /dev/null +++ b/pyomo/contrib/observer/component_collector.py @@ -0,0 +1,81 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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.expr.visitor import StreamBasedExpressionVisitor +from pyomo.core.expr.numeric_expr import ( + ExternalFunctionExpression, + NPV_ExternalFunctionExpression, +) +from pyomo.core.base.var import VarData, ScalarVar +from pyomo.core.base.param import ParamData, ScalarParam +from pyomo.core.base.expression import ExpressionData, ScalarExpression + + +def handle_var(node, collector): + collector.variables[id(node)] = node + return None + + +def handle_param(node, collector): + collector.params[id(node)] = node + return None + + +def handle_named_expression(node, collector): + collector.named_expressions[id(node)] = node + return None + + +def handle_external_function(node, collector): + collector.external_functions[id(node)] = node + return None + + +collector_handlers = { + VarData: handle_var, + ScalarVar: handle_var, + ParamData: handle_param, + ScalarParam: handle_param, + ExpressionData: handle_named_expression, + ScalarExpression: handle_named_expression, + ExternalFunctionExpression: handle_external_function, + NPV_ExternalFunctionExpression: handle_external_function, +} + + +class _ComponentFromExprCollector(StreamBasedExpressionVisitor): + def __init__(self, **kwds): + self.named_expressions = {} + self.variables = {} + self.params = {} + self.external_functions = {} + super().__init__(**kwds) + + def exitNode(self, node, data): + nt = type(node) + if nt in collector_handlers: + return collector_handlers[nt](node, self) + else: + return None + + +_visitor = _ComponentFromExprCollector() + + +def collect_components_from_expr(expr): + _visitor.__init__() + _visitor.walk_expression(expr) + return ( + list(_visitor.named_expressions.values()), + list(_visitor.variables.values()), + list(_visitor.params.values()), + list(_visitor.external_functions.values()), + ) diff --git a/pyomo/contrib/observer/model_observer.py b/pyomo/contrib/observer/model_observer.py new file mode 100644 index 00000000000..4ab52100376 --- /dev/null +++ b/pyomo/contrib/observer/model_observer.py @@ -0,0 +1,742 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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. +# __________________________________________________________________________ + +import abc +import datetime +from typing import List, Sequence, Optional + +from pyomo.common.config import ConfigDict, ConfigValue +from pyomo.core.base.constraint import ConstraintData, Constraint +from pyomo.core.base.sos import SOSConstraintData, SOSConstraint +from pyomo.core.base.var import VarData +from pyomo.core.base.param import ParamData, Param +from pyomo.core.base.objective import ObjectiveData +from pyomo.core.staleflag import StaleFlagManager +from pyomo.common.collections import ComponentMap +from pyomo.common.timing import HierarchicalTimer +from pyomo.contrib.solver.common.results import Results +from pyomo.contrib.solver.common.util import get_objective +from .component_collector import collect_components_from_expr +from pyomo.common.numeric_types import native_numeric_types + + +""" +The ModelChangeDetector is meant to be used to automatically identify changes +in a Pyomo model or block. Here is a list of changes that will be detected. +Note that inactive components (e.g., constraints) are treated as "removed". + - new constraints that have been added to the model + - constraints that have been removed from the model + - new variables that have been detected in new or modified constraints/objectives + - old variables that are no longer used in any constraints/objectives + - new parameters that have been detected in new or modified constraints/objectives + - old parameters that are no longer used in any constraints/objectives + - new objectives that have been added to the model + - objectives that have been removed from the model + - modified constraint expressions (relies on expressions being immutable) + - modified objective expressions (relies on expressions being immutable) + - modified objective sense + - changes to variable bounds, domains, and "fixed" flags + - changes to named expressions (relies on expressions being immutable) + - changes to parameter values and fixed variable values +""" + + +class AutoUpdateConfig(ConfigDict): + """ + Control which parts of the model are automatically checked and/or updated upon re-solve + """ + + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + if doc is None: + doc = 'Configuration options to detect changes in model between solves' + super().__init__( + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + + self.check_for_new_or_removed_constraints: bool = self.declare( + 'check_for_new_or_removed_constraints', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, new/old constraints will not be automatically detected on subsequent + solves. Use False only when manually updating the solver with opt.add_constraints() + and opt.remove_constraints() or when you are certain constraints are not being + added to/removed from the model.""", + ), + ) + self.check_for_new_objective: bool = self.declare( + 'check_for_new_objective', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, new/old objectives will not be automatically detected on subsequent + solves. Use False only when manually updating the solver with opt.set_objective() or + when you are certain objectives are not being added to / removed from the model.""", + ), + ) + self.update_constraints: bool = self.declare( + 'update_constraints', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, changes to existing constraints will not be automatically detected on + subsequent solves. This includes changes to the lower, body, and upper attributes of + constraints. Use False only when manually updating the solver with + opt.remove_constraints() and opt.add_constraints() or when you are certain constraints + are not being modified.""", + ), + ) + self.update_vars: bool = self.declare( + 'update_vars', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, changes to existing variables will not be automatically detected on + subsequent solves. This includes changes to the lb, ub, domain, and fixed + attributes of variables. Use False only when manually updating the observer with + opt.update_variables() or when you are certain variables are not being modified. + Note that changes to values of fixed variables is handled by + update_parameters_and_fixed_vars.""", + ), + ) + self.update_parameters: bool = self.declare( + 'update_parameters', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, changes to parameter values and fixed variable values will + not be automatically detected on subsequent solves. Use False only + when manually updating the observer with + opt.update_parameters_and_fixed_variables() or when you are certain + parameters are not being modified.""", + ), + ) + self.update_named_expressions: bool = self.declare( + 'update_named_expressions', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, changes to Expressions will not be automatically detected on + subsequent solves. Use False only when manually updating the solver with + opt.remove_constraints() and opt.add_constraints() or when you are certain + Expressions are not being modified.""", + ), + ) + self.update_objective: bool = self.declare( + 'update_objective', + ConfigValue( + domain=bool, + default=True, + description=""" + If False, changes to objectives will not be automatically detected on + subsequent solves. This includes the expr and sense attributes of objectives. Use + False only when manually updating the solver with opt.set_objective() or when you are + certain objectives are not being modified.""", + ), + ) + + +class Observer(abc.ABC): + @abc.abstractmethod + def add_variables(self, variables: List[VarData]): + pass + + @abc.abstractmethod + def add_parameters(self, params: List[ParamData]): + pass + + @abc.abstractmethod + def add_constraints(self, cons: List[ConstraintData]): + pass + + @abc.abstractmethod + def add_sos_constraints(self, cons: List[SOSConstraintData]): + pass + + @abc.abstractmethod + def set_objective(self, obj: Optional[ObjectiveData]): + pass + + @abc.abstractmethod + def remove_constraints(self, cons: List[ConstraintData]): + pass + + @abc.abstractmethod + def remove_sos_constraints(self, cons: List[SOSConstraintData]): + pass + + @abc.abstractmethod + def remove_variables(self, variables: List[VarData]): + pass + + @abc.abstractmethod + def remove_parameters(self, params: List[ParamData]): + pass + + @abc.abstractmethod + def update_variables(self, variables: List[VarData]): + pass + + @abc.abstractmethod + def update_parameters(self, params: List[ParamData]): + pass + + +class ModelChangeDetector: + def __init__(self, observers: Sequence[Observer], **kwds): + """ + Parameters + ---------- + observers: Sequence[Observer] + The objects to notify when changes are made to the model + """ + self._observers: List[Observer] = list(observers) + self._model = None + self._active_constraints = {} # maps constraint to expression + self._active_sos = {} + self._vars = {} # maps var id to (var, lb, ub, fixed, domain, value) + self._params = {} # maps param id to param + self._objective = None + self._objective_expr = None + self._objective_sense = None + self._named_expressions = ( + {} + ) # maps constraint to list of tuples (named_expr, named_expr.expr) + self._external_functions = ComponentMap() + self._obj_named_expressions = [] + self._referenced_variables = ( + {} + ) # var_id: [dict[constraints, None], dict[sos constraints, None], None or objective] + self._referenced_params = ( + {} + ) # param_id: [dict[constraints, None], dict[sos constraints, None], None or objective] + self._vars_referenced_by_con = {} + self._vars_referenced_by_obj = [] + self._params_referenced_by_con = {} + self._params_referenced_by_obj = [] + self._expr_types = None + self.config: AutoUpdateConfig = AutoUpdateConfig()( + value=kwds, preserve_implicit=True + ) + + def set_instance(self, model): + saved_config = self.config + self.__init__(observers=self._observers) + self.config = saved_config + self._model = model + self._add_block(model) + + def add_variables(self, variables: List[VarData]): + for v in variables: + if id(v) in self._referenced_variables: + raise ValueError(f'Variable {v.name} has already been added') + self._referenced_variables[id(v)] = [{}, {}, None] + self._vars[id(v)] = ( + v, + v._lb, + v._ub, + v.fixed, + v.domain.get_interval(), + v.value, + ) + for obs in self._observers: + obs.add_variables(variables) + + def add_parameters(self, params: List[ParamData]): + for p in params: + pid = id(p) + if pid in self._referenced_params: + raise ValueError(f'Parameter {p.name} has already been added') + self._referenced_params[pid] = [{}, {}, None] + self._params[id(p)] = (p, p.value) + for obs in self._observers: + obs.add_parameters(params) + + def _check_for_new_vars(self, variables: List[VarData]): + new_vars = {} + for v in variables: + v_id = id(v) + if v_id not in self._referenced_variables: + new_vars[v_id] = v + self.add_variables(list(new_vars.values())) + + def _check_to_remove_vars(self, variables: List[VarData]): + vars_to_remove = {} + for v in variables: + v_id = id(v) + ref_cons, ref_sos, ref_obj = self._referenced_variables[v_id] + if len(ref_cons) == 0 and len(ref_sos) == 0 and ref_obj is None: + vars_to_remove[v_id] = v + self.remove_variables(list(vars_to_remove.values())) + + def _check_for_new_params(self, params: List[ParamData]): + new_params = {} + for p in params: + pid = id(p) + if pid not in self._referenced_params: + new_params[pid] = p + self.add_parameters(list(new_params.values())) + + def _check_to_remove_params(self, params: List[ParamData]): + params_to_remove = {} + for p in params: + p_id = id(p) + ref_cons, ref_sos, ref_obj = self._referenced_params[p_id] + if len(ref_cons) == 0 and len(ref_sos) == 0 and ref_obj is None: + params_to_remove[p_id] = p + self.remove_parameters(list(params_to_remove.values())) + + def add_constraints(self, cons: List[ConstraintData]): + vars_to_check = [] + params_to_check = [] + for con in cons: + if con in self._active_constraints: + raise ValueError(f'Constraint {con.name} has already been added') + self._active_constraints[con] = con.expr + tmp = collect_components_from_expr(con.expr) + named_exprs, variables, parameters, external_functions = tmp + vars_to_check.extend(variables) + params_to_check.extend(parameters) + self._named_expressions[con] = [(e, e.expr) for e in named_exprs] + if len(external_functions) > 0: + self._external_functions[con] = external_functions + self._vars_referenced_by_con[con] = variables + self._params_referenced_by_con[con] = parameters + self._check_for_new_vars(vars_to_check) + self._check_for_new_params(params_to_check) + for con in cons: + variables = self._vars_referenced_by_con[con] + parameters = self._params_referenced_by_con[con] + for v in variables: + self._referenced_variables[id(v)][0][con] = None + for p in parameters: + self._referenced_params[id(p)][0][con] = None + for obs in self._observers: + obs.add_constraints(cons) + + def add_sos_constraints(self, cons: List[SOSConstraintData]): + vars_to_check = [] + params_to_check = [] + for con in cons: + if con in self._active_sos: + raise ValueError(f'Constraint {con.name} has already been added') + sos_items = list(con.get_items()) + self._active_sos[con] = ( + [i[0] for i in sos_items], + [i[1] for i in sos_items], + ) + variables = [] + params = [] + for v, p in sos_items: + variables.append(v) + if type(p) in native_numeric_types: + continue + if p.is_parameter_type(): + params.append(p) + vars_to_check.extend(variables) + params_to_check.extend(params) + self._named_expressions[con] = [] + self._vars_referenced_by_con[con] = variables + self._params_referenced_by_con[con] = params + self._check_for_new_vars(vars_to_check) + self._check_for_new_params(params_to_check) + for con in cons: + variables = self._vars_referenced_by_con[con] + params = self._params_referenced_by_con[con] + for v in variables: + self._referenced_variables[id(v)][1][con] = None + for p in params: + self._referenced_params[id(p)][1][con] = None + for obs in self._observers: + obs.add_sos_constraints(cons) + + def set_objective(self, obj: Optional[ObjectiveData]): + vars_to_remove_check = [] + params_to_remove_check = [] + if self._objective is not None: + for v in self._vars_referenced_by_obj: + self._referenced_variables[id(v)][2] = None + for p in self._params_referenced_by_obj: + self._referenced_params[id(p)][2] = None + vars_to_remove_check.extend(self._vars_referenced_by_obj) + params_to_remove_check.extend(self._params_referenced_by_obj) + self._external_functions.pop(self._objective, None) + if obj is not None: + self._objective = obj + self._objective_expr = obj.expr + self._objective_sense = obj.sense + tmp = collect_components_from_expr(obj.expr) + named_exprs, variables, parameters, external_functions = tmp + self._check_for_new_vars(variables) + self._check_for_new_params(parameters) + self._obj_named_expressions = [(i, i.expr) for i in named_exprs] + if len(external_functions) > 0: + self._external_functions[obj] = external_functions + self._vars_referenced_by_obj = variables + self._params_referenced_by_obj = parameters + for v in variables: + self._referenced_variables[id(v)][2] = obj + for p in parameters: + self._referenced_params[id(p)][2] = obj + else: + self._vars_referenced_by_obj = [] + self._params_referenced_by_obj = [] + self._objective = None + self._objective_expr = None + self._objective_sense = None + self._obj_named_expressions = [] + for obs in self._observers: + obs.set_objective(obj) + self._check_to_remove_vars(vars_to_remove_check) + self._check_to_remove_params(params_to_remove_check) + + def _add_block(self, block): + self.add_constraints( + list( + block.component_data_objects(Constraint, descend_into=True, active=True) + ) + ) + self.add_sos_constraints( + list( + block.component_data_objects( + SOSConstraint, descend_into=True, active=True + ) + ) + ) + obj = get_objective(block) + self.set_objective(obj) + + def remove_constraints(self, cons: List[ConstraintData]): + for obs in self._observers: + obs.remove_constraints(cons) + vars_to_check = [] + params_to_check = [] + for con in cons: + if con not in self._active_constraints: + raise ValueError( + f'Cannot remove constraint {con.name} - it was not added' + ) + for v in self._vars_referenced_by_con[con]: + self._referenced_variables[id(v)][0].pop(con) + for p in self._params_referenced_by_con[con]: + self._referenced_params[id(p)][0].pop(con) + vars_to_check.extend(self._vars_referenced_by_con[con]) + params_to_check.extend(self._params_referenced_by_con[con]) + del self._active_constraints[con] + del self._named_expressions[con] + self._external_functions.pop(con, None) + del self._vars_referenced_by_con[con] + del self._params_referenced_by_con[con] + self._check_to_remove_vars(vars_to_check) + self._check_to_remove_params(params_to_check) + + def remove_sos_constraints(self, cons: List[SOSConstraintData]): + for obs in self._observers: + obs.remove_sos_constraints(cons) + vars_to_check = [] + params_to_check = [] + for con in cons: + if con not in self._active_sos: + raise ValueError( + f'Cannot remove constraint {con.name} - it was not added' + ) + for v in self._vars_referenced_by_con[con]: + self._referenced_variables[id(v)][1].pop(con) + for p in self._params_referenced_by_con[con]: + self._referenced_params[id(p)][1].pop(con) + vars_to_check.extend(self._vars_referenced_by_con[con]) + params_to_check.extend(self._params_referenced_by_con[con]) + del self._active_sos[con] + del self._named_expressions[con] + del self._vars_referenced_by_con[con] + del self._params_referenced_by_con[con] + self._check_to_remove_vars(vars_to_check) + self._check_to_remove_params(params_to_check) + + def remove_variables(self, variables: List[VarData]): + for obs in self._observers: + obs.remove_variables(variables) + for v in variables: + v_id = id(v) + if v_id not in self._referenced_variables: + raise ValueError( + f'Cannot remove variable {v.name} - it has not been added' + ) + cons_using, sos_using, obj_using = self._referenced_variables[v_id] + if cons_using or sos_using or (obj_using is not None): + raise ValueError( + f'Cannot remove variable {v.name} - it is still being used by constraints or the objective' + ) + del self._referenced_variables[v_id] + del self._vars[v_id] + + def remove_parameters(self, params: List[ParamData]): + for obs in self._observers: + obs.remove_parameters(params) + for p in params: + p_id = id(p) + if p_id not in self._referenced_params: + raise ValueError( + f'Cannot remove parameter {p.name} - it has not been added' + ) + cons_using, sos_using, obj_using = self._referenced_params[p_id] + if cons_using or sos_using or (obj_using is not None): + raise ValueError( + f'Cannot remove parameter {p.name} - it is still being used by constraints or the objective' + ) + del self._referenced_params[p_id] + del self._params[p_id] + + def update_variables(self, variables: List[VarData]): + for v in variables: + self._vars[id(v)] = ( + v, + v._lb, + v._ub, + v.fixed, + v.domain.get_interval(), + v.value, + ) + for obs in self._observers: + obs.update_variables(variables) + + def update_parameters(self, params): + for p in params: + self._params[id(p)] = (p, p.value) + for obs in self._observers: + obs.update_parameters(params) + + def _check_for_new_or_removed_sos(self): + new_sos = [] + old_sos = [] + current_sos_dict = { + c: None + for c in self._model.component_data_objects( + SOSConstraint, descend_into=True, active=True + ) + } + for c in current_sos_dict.keys(): + if c not in self._active_sos: + new_sos.append(c) + for c in self._active_sos: + if c not in current_sos_dict: + old_sos.append(c) + return new_sos, old_sos + + def _check_for_new_or_removed_constraints(self): + new_cons = [] + old_cons = [] + current_cons_dict = { + c: None + for c in self._model.component_data_objects( + Constraint, descend_into=True, active=True + ) + } + for c in current_cons_dict.keys(): + if c not in self._active_constraints: + new_cons.append(c) + for c in self._active_constraints: + if c not in current_cons_dict: + old_cons.append(c) + return new_cons, old_cons + + def _check_for_modified_sos(self): + sos_to_update = [] + for c, (old_vlist, old_plist) in self._active_sos.items(): + sos_items = list(c.get_items()) + new_vlist = [i[0] for i in sos_items] + new_plist = [i[1] for i in sos_items] + if len(old_vlist) != len(new_vlist): + sos_to_update.append(c) + elif len(old_plist) != len(new_plist): + sos_to_update.append(c) + else: + needs_update = False + for v1, v2 in zip(old_vlist, new_vlist): + if v1 is not v2: + needs_update = True + break + for p1, p2 in zip(old_plist, new_plist): + if p1 is not p2: + needs_update = True + if needs_update: + break + if needs_update: + sos_to_update.append(c) + return sos_to_update + + def _check_for_modified_constraints(self): + cons_to_update = [] + for c, expr in self._active_constraints.items(): + if c.expr is not expr: + cons_to_update.append(c) + return cons_to_update + + def _check_for_var_changes(self): + vars_to_update = [] + cons_to_update = {} + update_obj = False + for vid, (v, _lb, _ub, _fixed, _domain_interval, _value) in self._vars.items(): + if v.fixed != _fixed: + vars_to_update.append(v) + for c in self._referenced_variables[vid][0]: + cons_to_update[c] = None + if self._referenced_variables[vid][2] is not None: + update_obj = True + elif v._lb is not _lb: + vars_to_update.append(v) + elif v._ub is not _ub: + vars_to_update.append(v) + elif _domain_interval != v.domain.get_interval(): + vars_to_update.append(v) + elif v.value != _value: + vars_to_update.append(v) + cons_to_update = list(cons_to_update.keys()) + return vars_to_update, cons_to_update, update_obj + + def _check_for_param_changes(self): + params_to_update = [] + for pid, (p, val) in self._params.items(): + if p.value != val: + params_to_update.append(p) + return params_to_update + + def _check_for_named_expression_changes(self): + cons_to_update = [] + for con, ne_list in self._named_expressions.items(): + for named_expr, old_expr in ne_list: + if named_expr.expr is not old_expr: + cons_to_update.append(con) + break + update_obj = False + ne_list = self._obj_named_expressions + for named_expr, old_expr in ne_list: + if named_expr.expr is not old_expr: + update_obj = True + break + return cons_to_update, update_obj + + def _check_for_new_objective(self): + update_obj = False + new_obj = get_objective(self._model) + if new_obj is not self._objective: + update_obj = True + return new_obj, update_obj + + def _check_for_objective_changes(self): + update_obj = False + if self._objective is None: + return update_obj + if self._objective.expr is not self._objective_expr: + update_obj = True + elif self._objective.sense != self._objective_sense: + # we can definitely do something faster here than resetting the whole objective + update_obj = True + return update_obj + + def update(self, timer: Optional[HierarchicalTimer] = None, **kwds): + if timer is None: + timer = HierarchicalTimer() + config: AutoUpdateConfig = self.config(value=kwds, preserve_implicit=True) + + added_cons = set() + added_sos = set() + + if config.check_for_new_or_removed_constraints: + timer.start('sos') + new_sos, old_sos = self._check_for_new_or_removed_sos() + self.add_sos_constraints(new_sos) + self.remove_sos_constraints(old_sos) + added_sos.update(new_sos) + timer.stop('sos') + timer.start('cons') + new_cons, old_cons = self._check_for_new_or_removed_constraints() + self.add_constraints(new_cons) + self.remove_constraints(old_cons) + added_cons.update(new_cons) + timer.stop('cons') + + if config.update_constraints: + timer.start('cons') + cons_to_update = self._check_for_modified_constraints() + self.remove_constraints(cons_to_update) + self.add_constraints(cons_to_update) + added_cons.update(cons_to_update) + timer.stop('cons') + timer.start('sos') + sos_to_update = self._check_for_modified_sos() + self.remove_sos_constraints(sos_to_update) + self.add_sos_constraints(sos_to_update) + added_sos.update(sos_to_update) + timer.stop('sos') + + need_to_set_objective = False + + if config.update_vars: + timer.start('vars') + vars_to_update, cons_to_update, update_obj = self._check_for_var_changes() + self.update_variables(vars_to_update) + cons_to_update = [i for i in cons_to_update if i not in added_cons] + self.remove_constraints(cons_to_update) + self.add_constraints(cons_to_update) + added_cons.update(cons_to_update) + if update_obj: + need_to_set_objective = True + timer.stop('vars') + + if config.update_named_expressions: + timer.start('named expressions') + cons_to_update, update_obj = self._check_for_named_expression_changes() + cons_to_update = [i for i in cons_to_update if i not in added_cons] + self.remove_constraints(cons_to_update) + self.add_constraints(cons_to_update) + added_cons.update(cons_to_update) + if update_obj: + need_to_set_objective = True + timer.stop('named expressions') + + timer.start('objective') + new_obj = self._objective + if config.check_for_new_objective: + new_obj, update_obj = self._check_for_new_objective() + if update_obj: + need_to_set_objective = True + if config.update_objective: + update_obj = self._check_for_objective_changes() + if update_obj: + need_to_set_objective = True + + if need_to_set_objective: + self.set_objective(new_obj) + timer.stop('objective') + + if config.update_parameters: + timer.start('params') + params_to_update = self._check_for_param_changes() + self.update_parameters(params_to_update) + timer.stop('params') diff --git a/pyomo/contrib/observer/tests/__init__.py b/pyomo/contrib/observer/tests/__init__.py new file mode 100644 index 00000000000..6eb9ea8b81d --- /dev/null +++ b/pyomo/contrib/observer/tests/__init__.py @@ -0,0 +1,10 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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. +# ___________________________________________________________________________ diff --git a/pyomo/contrib/observer/tests/test_change_detector.py b/pyomo/contrib/observer/tests/test_change_detector.py new file mode 100644 index 00000000000..f3ce6e28b5c --- /dev/null +++ b/pyomo/contrib/observer/tests/test_change_detector.py @@ -0,0 +1,365 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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.constraint import ConstraintData +from pyomo.core.base.objective import ObjectiveData +from pyomo.core.base.param import ParamData +from pyomo.core.base.sos import SOSConstraintData +from pyomo.core.base.var import VarData +import pyomo.environ as pyo +from pyomo.common import unittest +from typing import List +from pyomo.contrib.observer.model_observer import ( + Observer, + ModelChangeDetector, + AutoUpdateConfig, +) +from pyomo.common.collections import ComponentMap +import logging + + +logger = logging.getLogger(__name__) + + +def make_count_dict(): + d = {'add': 0, 'remove': 0, 'update': 0, 'set': 0} + return d + + +class ObserverChecker(Observer): + def __init__(self): + super().__init__() + self.counts = ComponentMap() + """ + counts is a mapping from component (e.g., variable) to another + mapping from string ('add', 'remove', 'update', or 'set') to an int that + indicates the number of times the corresponding method has been called + """ + + def check(self, expected): + unittest.assertStructuredAlmostEqual( + first=expected, second=self.counts, places=7 + ) + + def _process(self, comps, key): + for c in comps: + if c not in self.counts: + self.counts[c] = make_count_dict() + self.counts[c][key] += 1 + + def pprint(self): + for k, d in self.counts.items(): + print(f'{k}:') + for a, v in d.items(): + print(f' {a}: {v}') + + def add_variables(self, variables: List[VarData]): + for v in variables: + assert v.is_variable_type() + self._process(variables, 'add') + + def add_parameters(self, params: List[ParamData]): + for p in params: + assert p.is_parameter_type() + self._process(params, 'add') + + def add_constraints(self, cons: List[ConstraintData]): + for c in cons: + assert isinstance(c, ConstraintData) + self._process(cons, 'add') + + def add_sos_constraints(self, cons: List[SOSConstraintData]): + for c in cons: + assert isinstance(c, SOSConstraintData) + self._process(cons, 'add') + + def set_objective(self, obj: ObjectiveData): + assert obj is None or isinstance(obj, ObjectiveData) + self._process([obj], 'set') + + def remove_constraints(self, cons: List[ConstraintData]): + for c in cons: + assert isinstance(c, ConstraintData) + self._process(cons, 'remove') + + def remove_sos_constraints(self, cons: List[SOSConstraintData]): + for c in cons: + assert isinstance(c, SOSConstraintData) + self._process(cons, 'remove') + + def remove_variables(self, variables: List[VarData]): + for v in variables: + assert v.is_variable_type() + self._process(variables, 'remove') + + def remove_parameters(self, params: List[ParamData]): + for p in params: + assert p.is_parameter_type() + self._process(params, 'remove') + + def update_variables(self, variables: List[VarData]): + for v in variables: + assert v.is_variable_type() + self._process(variables, 'update') + + def update_parameters(self, params: List[ParamData]): + for p in params: + assert p.is_parameter_type() + self._process(params, 'update') + + +class TestChangeDetector(unittest.TestCase): + def test_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.p = pyo.Param(mutable=True, initialize=1) + + obs = ObserverChecker() + detector = ModelChangeDetector([obs]) + + expected = ComponentMap() + expected[None] = make_count_dict() + expected[None]['set'] += 1 + + detector.set_instance(m) + obs.check(expected) + + m.obj = pyo.Objective(expr=m.x**2 + m.p * m.y**2) + detector.update() + expected[m.obj] = make_count_dict() + expected[m.obj]['set'] += 1 + expected[m.x] = make_count_dict() + expected[m.x]['add'] += 1 + expected[m.y] = make_count_dict() + expected[m.y]['add'] += 1 + expected[m.p] = make_count_dict() + expected[m.p]['add'] += 1 + obs.check(expected) + + m.y.setlb(0) + detector.update() + expected[m.y]['update'] += 1 + obs.check(expected) + + m.x.fix(2) + detector.update() + expected[m.x]['update'] += 1 + expected[m.obj]['set'] += 1 + obs.check(expected) + + m.x.unfix() + detector.update() + expected[m.x]['update'] += 1 + expected[m.obj]['set'] += 1 + obs.check(expected) + + m.p.value = 2 + detector.update() + expected[m.p]['update'] += 1 + obs.check(expected) + + m.obj.expr = m.x**2 + m.y**2 + detector.update() + expected[m.p]['remove'] += 1 + expected[m.obj]['set'] += 1 + obs.check(expected) + + del m.obj + m.obj = pyo.Objective(expr=m.p * m.x) + detector.update() + expected[m.p]['add'] += 1 + expected[m.y]['remove'] += 1 + # remember, m.obj is a different object now + expected[m.obj] = make_count_dict() + expected[m.obj]['set'] += 1 + + def test_constraints(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.p = pyo.Param(mutable=True, initialize=1) + + obs = ObserverChecker() + detector = ModelChangeDetector([obs]) + + expected = ComponentMap() + expected[None] = make_count_dict() + expected[None]['set'] += 1 + + detector.set_instance(m) + obs.check(expected) + + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= (m.x - m.p) ** 2) + detector.update() + expected[m.x] = make_count_dict() + expected[m.y] = make_count_dict() + expected[m.p] = make_count_dict() + expected[m.x]['add'] += 1 + expected[m.y]['add'] += 1 + expected[m.p]['add'] += 1 + expected[m.c1] = make_count_dict() + expected[m.c1]['add'] += 1 + expected[m.obj] = make_count_dict() + expected[m.obj]['set'] += 1 + obs.check(expected) + + # now fix a variable and make sure the + # constraint gets removed and added + m.x.fix(1) + obs.pprint() + detector.update() + obs.pprint() + expected[m.c1]['remove'] += 1 + expected[m.c1]['add'] += 1 + # because x and p are only used in the + # one constraint, they get removed when + # the constraint is removed and then + # added again when the constraint is added + expected[m.x]['update'] += 1 + expected[m.x]['remove'] += 1 + expected[m.x]['add'] += 1 + expected[m.p]['remove'] += 1 + expected[m.p]['add'] += 1 + obs.check(expected) + + def test_sos(self): + m = pyo.ConcreteModel() + m.a = pyo.Set(initialize=[1, 2, 3], ordered=True) + m.x = pyo.Var(m.a, within=pyo.Binary) + m.y = pyo.Var(within=pyo.Binary) + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.SOSConstraint(var=m.x, sos=1) + + obs = ObserverChecker() + detector = ModelChangeDetector([obs]) + detector.set_instance(m) + + expected = ComponentMap() + expected[m.obj] = make_count_dict() + for i in m.a: + expected[m.x[i]] = make_count_dict() + expected[m.y] = make_count_dict() + expected[m.c1] = make_count_dict() + expected[m.obj]['set'] += 1 + for i in m.a: + expected[m.x[i]]['add'] += 1 + expected[m.y]['add'] += 1 + expected[m.c1]['add'] += 1 + obs.check(expected) + + for i in m.a: + expected[m.x[i]]['remove'] += 1 + expected[m.c1]['remove'] += 1 + del m.c1 + detector.update() + obs.check(expected) + + def test_vars_and_params_elsewhere(self): + m1 = pyo.ConcreteModel() + m1.x = pyo.Var() + m1.y = pyo.Var() + m1.p = pyo.Param(mutable=True, initialize=1) + + m2 = pyo.ConcreteModel() + + obs = ObserverChecker() + detector = ModelChangeDetector([obs]) + + expected = ComponentMap() + expected[None] = make_count_dict() + expected[None]['set'] += 1 + + detector.set_instance(m2) + obs.check(expected) + + m2.obj = pyo.Objective(expr=m1.y) + m2.c1 = pyo.Constraint(expr=m1.y >= (m1.x - m1.p) ** 2) + detector.update() + expected[m1.x] = make_count_dict() + expected[m1.y] = make_count_dict() + expected[m1.p] = make_count_dict() + expected[m1.x]['add'] += 1 + expected[m1.y]['add'] += 1 + expected[m1.p]['add'] += 1 + expected[m2.c1] = make_count_dict() + expected[m2.c1]['add'] += 1 + expected[m2.obj] = make_count_dict() + expected[m2.obj]['set'] += 1 + obs.check(expected) + + # now fix a variable and make sure the + # constraint gets removed and added + m1.x.fix(1) + obs.pprint() + detector.update() + obs.pprint() + expected[m2.c1]['remove'] += 1 + expected[m2.c1]['add'] += 1 + # because x and p are only used in the + # one constraint, they get removed when + # the constraint is removed and then + # added again when the constraint is added + expected[m1.x]['update'] += 1 + expected[m1.x]['remove'] += 1 + expected[m1.x]['add'] += 1 + expected[m1.p]['remove'] += 1 + expected[m1.p]['add'] += 1 + obs.check(expected) + + def test_named_expression(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.p = pyo.Param(mutable=True, initialize=1) + + obs = ObserverChecker() + detector = ModelChangeDetector([obs]) + + expected = ComponentMap() + expected[None] = make_count_dict() + expected[None]['set'] += 1 + + detector.set_instance(m) + obs.check(expected) + + m.obj = pyo.Objective(expr=m.y) + m.e = pyo.Expression(expr=m.x - m.p) + m.c1 = pyo.Constraint(expr=m.y >= m.e) + detector.update() + expected[m.x] = make_count_dict() + expected[m.y] = make_count_dict() + expected[m.p] = make_count_dict() + expected[m.x]['add'] += 1 + expected[m.y]['add'] += 1 + expected[m.p]['add'] += 1 + expected[m.c1] = make_count_dict() + expected[m.c1]['add'] += 1 + expected[m.obj] = make_count_dict() + expected[m.obj]['set'] += 1 + obs.check(expected) + + # now modify the named expression and make sure the + # constraint gets removed and added + m.e.expr = (m.x - m.p) ** 2 + detector.update() + expected[m.c1]['remove'] += 1 + expected[m.c1]['add'] += 1 + # because x and p are only used in the + # one constraint, they get removed when + # the constraint is removed and then + # added again when the constraint is added + expected[m.x]['remove'] += 1 + expected[m.x]['add'] += 1 + expected[m.p]['remove'] += 1 + expected[m.p]['add'] += 1 + obs.check(expected) diff --git a/pyomo/contrib/solver/plugins.py b/pyomo/contrib/solver/plugins.py index 19fc9b2b2a1..bba5dc4db3f 100644 --- a/pyomo/contrib/solver/plugins.py +++ b/pyomo/contrib/solver/plugins.py @@ -12,8 +12,8 @@ from .common.factory import SolverFactory from .solvers.ipopt import Ipopt, LegacyIpoptSolver -from .solvers.gurobi_persistent import GurobiPersistent -from .solvers.gurobi_direct import GurobiDirect +from .solvers.gurobi.gurobi_direct import GurobiDirect +from .solvers.gurobi.gurobi_persistent import GurobiDirectQuadratic, GurobiPersistent from .solvers.highs import Highs @@ -32,5 +32,10 @@ def load(): doc='Direct (scipy-based) interface to Gurobi', )(GurobiDirect) SolverFactory.register( - name='highs', legacy_name='highs', doc='Persistent interface to HiGHS' + name='gurobi_direct_quadratic', + legacy_name='gurobi_direct_quadratic_v2', + doc='Direct interface to Gurobi', + )(GurobiDirectQuadratic) + SolverFactory.register( + name='highs', legacy_name='highs_v2', doc='Persistent interface to HiGHS' )(Highs) diff --git a/pyomo/contrib/solver/solvers/gurobi/__init__.py b/pyomo/contrib/solver/solvers/gurobi/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pyomo/contrib/solver/solvers/gurobi/gurobi_direct.py b/pyomo/contrib/solver/solvers/gurobi/gurobi_direct.py new file mode 100644 index 00000000000..16c633c7d7c --- /dev/null +++ b/pyomo/contrib/solver/solvers/gurobi/gurobi_direct.py @@ -0,0 +1,184 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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. +# ___________________________________________________________________________ + +import operator + +from pyomo.common.collections import ComponentMap, ComponentSet +from pyomo.common.shutdown import python_is_shutting_down +from pyomo.core.staleflag import StaleFlagManager +from pyomo.repn.plugins.standard_form import LinearStandardFormCompiler + +from pyomo.contrib.solver.common.util import ( + NoDualsError, + NoReducedCostsError, + NoSolutionError, + IncompatibleModelError, +) +from pyomo.contrib.solver.common.solution_loader import SolutionLoaderBase +from .gurobi_direct_base import GurobiDirectBase, gurobipy + + +class GurobiDirectSolutionLoader(SolutionLoaderBase): + def __init__(self, grb_model, grb_cons, grb_vars, pyo_cons, pyo_vars): + self._grb_model = grb_model + self._grb_cons = grb_cons + self._grb_vars = grb_vars + self._pyo_cons = pyo_cons + self._pyo_vars = pyo_vars + GurobiDirectBase._register_env_client() + + def __del__(self): + if python_is_shutting_down(): + return + # Free the associated model + if self._grb_model is not None: + self._grb_cons = None + self._grb_vars = None + self._pyo_cons = None + self._pyo_vars = None + # explicitly release the model + self._grb_model.dispose() + self._grb_model = None + # Release the gurobi license if this is the last reference to + # the environment (either through a results object or solver + # interface) + GurobiDirectBase._release_env_client() + + def load_vars(self, vars_to_load=None, solution_number=0): + assert solution_number == 0 + if self._grb_model.SolCount == 0: + raise NoSolutionError() + + iterator = zip(self._pyo_vars, self._grb_vars.x.tolist()) + if vars_to_load: + vars_to_load = ComponentSet(vars_to_load) + iterator = filter(lambda var_val: var_val[0] in vars_to_load, iterator) + for p_var, g_var in iterator: + p_var.set_value(g_var, skip_validation=True) + StaleFlagManager.mark_all_as_stale(delayed=True) + + def get_primals(self, vars_to_load=None, solution_number=0): + assert solution_number == 0 + if self._grb_model.SolCount == 0: + raise NoSolutionError() + + iterator = zip(self._pyo_vars, self._grb_vars.x.tolist()) + if vars_to_load: + vars_to_load = ComponentSet(vars_to_load) + iterator = filter(lambda var_val: var_val[0] in vars_to_load, iterator) + return ComponentMap(iterator) + + def get_duals(self, cons_to_load=None): + if self._grb_model.Status != gurobipy.GRB.OPTIMAL: + raise NoDualsError() + + def dedup(_iter): + last = None + for con_info_dual in _iter: + if not con_info_dual[1] and con_info_dual[0][0] is last: + continue + last = con_info_dual[0][0] + yield con_info_dual + + iterator = dedup(zip(self._pyo_cons, self._grb_cons.getAttr('Pi').tolist())) + if cons_to_load: + cons_to_load = set(cons_to_load) + iterator = filter( + lambda con_info_dual: con_info_dual[0][0] in cons_to_load, iterator + ) + return {con_info[0]: dual for con_info, dual in iterator} + + def get_reduced_costs(self, vars_to_load=None): + if self._grb_model.Status != gurobipy.GRB.OPTIMAL: + raise NoReducedCostsError() + + iterator = zip(self._pyo_vars, self._grb_vars.getAttr('Rc').tolist()) + if vars_to_load: + vars_to_load = ComponentSet(vars_to_load) + iterator = filter(lambda var_rc: var_rc[0] in vars_to_load, iterator) + return ComponentMap(iterator) + + +class GurobiDirect(GurobiDirectBase): + _minimum_version = (9, 0, 0) + + def __init__(self, **kwds): + super().__init__(**kwds) + self._gurobi_vars = None + self._pyomo_vars = None + + def _pyomo_gurobi_var_iter(self): + return zip(self._pyomo_vars, self._gurobi_vars.tolist()) + + def _create_solver_model(self, pyomo_model): + timer = self.config.timer + + timer.start('compile_model') + repn = LinearStandardFormCompiler().write( + pyomo_model, mixed_form=True, set_sense=None + ) + timer.stop('compile_model') + + if len(repn.objectives) > 1: + raise IncompatibleModelError( + f"The {self.__class__.__name__} solver only supports models " + f"with zero or one objectives (received {len(repn.objectives)})." + ) + + timer.start('prepare_matrices') + inf = float('inf') + ninf = -inf + bounds = list(map(operator.attrgetter('bounds'), repn.columns)) + lb = [ninf if _b is None else _b for _b in map(operator.itemgetter(0), bounds)] + ub = [inf if _b is None else _b for _b in map(operator.itemgetter(1), bounds)] + CON = gurobipy.GRB.CONTINUOUS + BIN = gurobipy.GRB.BINARY + INT = gurobipy.GRB.INTEGER + vtype = [ + ( + CON + if v.is_continuous() + else BIN if v.is_binary() else INT if v.is_integer() else '?' + ) + for v in repn.columns + ] + sense_type = list('=<>') # Note: ordering matches 0, 1, -1 + sense = [sense_type[r[1]] for r in repn.rows] + timer.stop('prepare_matrices') + + gurobi_model = gurobipy.Model(env=self.env()) + + timer.start('transfer_model') + x = gurobi_model.addMVar( + len(repn.columns), + lb=lb, + ub=ub, + obj=repn.c.todense()[0] if repn.c.shape[0] else 0, + vtype=vtype, + ) + A = gurobi_model.addMConstr(repn.A, x, sense, repn.rhs) + if repn.c.shape[0]: + gurobi_model.setAttr('ObjCon', repn.c_offset[0]) + gurobi_model.setAttr('ModelSense', int(repn.objectives[0].sense)) + # Note: calling gurobi_model.update() here is not + # necessary (it will happen as part of optimize()): + # gurobi_model.update() + timer.stop('transfer_model') + + self._pyomo_vars = repn.columns + self._gurobi_vars = x + + solution_loader = GurobiDirectSolutionLoader( + gurobi_model, A, x, repn.rows, repn.columns + ) + has_obj = len(repn.objectives) > 0 + + return gurobi_model, solution_loader, has_obj diff --git a/pyomo/contrib/solver/solvers/gurobi_direct.py b/pyomo/contrib/solver/solvers/gurobi/gurobi_direct_base.py similarity index 50% rename from pyomo/contrib/solver/solvers/gurobi_direct.py rename to pyomo/contrib/solver/solvers/gurobi/gurobi_direct_base.py index 45ea9dcc873..df6bb8b5327 100644 --- a/pyomo/contrib/solver/solvers/gurobi_direct.py +++ b/pyomo/contrib/solver/solvers/gurobi/gurobi_direct_base.py @@ -12,19 +12,17 @@ import datetime import io import math -import operator import os -from pyomo.common.collections import ComponentMap, ComponentSet +from pyomo.common.collections import ComponentMap from pyomo.common.config import ConfigValue from pyomo.common.dependencies import attempt_import from pyomo.common.enums import ObjectiveSense -from pyomo.common.errors import MouseTrap, ApplicationError +from pyomo.common.errors import ApplicationError from pyomo.common.shutdown import python_is_shutting_down from pyomo.common.tee import capture_output, TeeStream from pyomo.common.timing import HierarchicalTimer from pyomo.core.staleflag import StaleFlagManager -from pyomo.repn.plugins.standard_form import LinearStandardFormCompiler from pyomo.contrib.solver.common.base import SolverBase, Availability from pyomo.contrib.solver.common.config import BranchAndBoundConfig @@ -34,37 +32,22 @@ NoDualsError, NoReducedCostsError, NoSolutionError, - IncompatibleModelError, ) from pyomo.contrib.solver.common.results import ( Results, SolutionStatus, TerminationCondition, ) -from pyomo.contrib.solver.common.solution_loader import SolutionLoaderBase +import logging -gurobipy, gurobipy_available = attempt_import('gurobipy') - +logger = logging.getLogger(__name__) -class GurobiConfigMixin: - """ - Mixin class for Gurobi-specific configurations - """ - def __init__(self): - self.use_mipstart: bool = self.declare( - 'use_mipstart', - ConfigValue( - default=False, - domain=bool, - description="If True, the current values of the integer variables " - "will be passed to Gurobi.", - ), - ) +gurobipy, gurobipy_available = attempt_import('gurobipy') -class GurobiConfig(BranchAndBoundConfig, GurobiConfigMixin): +class GurobiConfig(BranchAndBoundConfig): def __init__( self, description=None, @@ -81,103 +64,150 @@ def __init__( implicit_domain=implicit_domain, visibility=visibility, ) - GurobiConfigMixin.__init__(self) + self.use_mipstart: bool = self.declare( + 'use_mipstart', + ConfigValue( + default=False, + domain=bool, + description="If True, the current values of the integer variables " + "will be passed to Gurobi.", + ), + ) -class GurobiDirectSolutionLoader(SolutionLoaderBase): - def __init__(self, grb_model, grb_cons, grb_vars, pyo_cons, pyo_vars, pyo_obj): - self._grb_model = grb_model - self._grb_cons = grb_cons - self._grb_vars = grb_vars - self._pyo_cons = pyo_cons - self._pyo_vars = pyo_vars - self._pyo_obj = pyo_obj - GurobiDirect._register_env_client() +def _load_suboptimal_mip_solution(solver_model, var_map, vars_to_load, solution_number): + """ + solver_model: gurobipy.Model + var_map: Dict[int, gurobipy.Var] + Maps the id of the pyomo variable to the gurobipy variable + vars_to_load: List[VarData] + solution_number: int + """ + if ( + solver_model.getAttr('NumIntVars') == 0 + and solver_model.getAttr('NumBinVars') == 0 + ): + raise ValueError('Cannot obtain suboptimal solutions for a continuous model') + original_solution_number = solver_model.getParamInfo('SolutionNumber')[2] + solver_model.setParam('SolutionNumber', solution_number) + gurobi_vars_to_load = [var_map[id(v)] for v in vars_to_load] + vals = solver_model.getAttr("Xn", gurobi_vars_to_load) + res = ComponentMap() + for var, val in zip(vars_to_load, vals): + res[var] = val + solver_model.setParam('SolutionNumber', original_solution_number) + return res + + +def _load_vars(solver_model, var_map, vars_to_load, solution_number=0): + """ + solver_model: gurobipy.Model + var_map: Dict[int, gurobipy.Var] + Maps the id of the pyomo variable to the gurobipy variable + vars_to_load: List[VarData] + solution_number: int + """ + for v, val in _get_primals( + solver_model=solver_model, + var_map=var_map, + vars_to_load=vars_to_load, + solution_number=solution_number, + ).items(): + v.set_value(val, skip_validation=True) + StaleFlagManager.mark_all_as_stale(delayed=True) - def __del__(self): - if python_is_shutting_down(): - return - # Free the associated model - if self._grb_model is not None: - self._grb_cons = None - self._grb_vars = None - self._pyo_cons = None - self._pyo_vars = None - self._pyo_obj = None - # explicitly release the model - self._grb_model.dispose() - self._grb_model = None - # Release the gurobi license if this is the last reference to - # the environment (either through a results object or solver - # interface) - GurobiDirect._release_env_client() - - def load_vars(self, vars_to_load=None, solution_number=0): - assert solution_number == 0 - if self._grb_model.SolCount == 0: - raise NoSolutionError() - - iterator = zip(self._pyo_vars, self._grb_vars.x.tolist()) - if vars_to_load: - vars_to_load = ComponentSet(vars_to_load) - iterator = filter(lambda var_val: var_val[0] in vars_to_load, iterator) - for p_var, g_var in iterator: - p_var.set_value(g_var, skip_validation=True) - StaleFlagManager.mark_all_as_stale(delayed=True) - - def get_primals(self, vars_to_load=None, solution_number=0): - assert solution_number == 0 - if self._grb_model.SolCount == 0: - raise NoSolutionError() - - iterator = zip(self._pyo_vars, self._grb_vars.x.tolist()) - if vars_to_load: - vars_to_load = ComponentSet(vars_to_load) - iterator = filter(lambda var_val: var_val[0] in vars_to_load, iterator) - return ComponentMap(iterator) - - def get_duals(self, cons_to_load=None): - if self._grb_model.Status != gurobipy.GRB.OPTIMAL: - raise NoDualsError() - - def dedup(_iter): - last = None - for con_info_dual in _iter: - if not con_info_dual[1] and con_info_dual[0][0] is last: - continue - last = con_info_dual[0][0] - yield con_info_dual - - iterator = dedup(zip(self._pyo_cons, self._grb_cons.getAttr('Pi').tolist())) - if cons_to_load: - cons_to_load = set(cons_to_load) - iterator = filter( - lambda con_info_dual: con_info_dual[0][0] in cons_to_load, iterator - ) - return {con_info[0]: dual for con_info, dual in iterator} - def get_reduced_costs(self, vars_to_load=None): - if self._grb_model.Status != gurobipy.GRB.OPTIMAL: - raise NoReducedCostsError() +def _get_primals(solver_model, var_map, vars_to_load, solution_number=0): + """ + solver_model: gurobipy.Model + var_map: Dict[int, gurobipy.Var] + Maps the id of the pyomo variable to the gurobipy variable + vars_to_load: List[VarData] + solution_number: int + """ + if solver_model.SolCount == 0: + raise NoSolutionError() + + if solution_number != 0: + return _load_suboptimal_mip_solution( + solver_model=solver_model, + var_map=var_map, + vars_to_load=vars_to_load, + solution_number=solution_number, + ) + + gurobi_vars_to_load = [var_map[id(v)] for v in vars_to_load] + vals = solver_model.getAttr("X", gurobi_vars_to_load) + + res = ComponentMap() + for var, val in zip(vars_to_load, vals): + res[var] = val + return res + + +def _get_reduced_costs(solver_model, var_map, vars_to_load): + """ + solver_model: gurobipy.Model + var_map: Dict[int, gurobipy.Var] + Maps the id of the pyomo variable to the gurobipy variable + vars_to_load: List[VarData] + """ + if solver_model.Status != gurobipy.GRB.OPTIMAL: + raise NoReducedCostsError() - iterator = zip(self._pyo_vars, self._grb_vars.getAttr('Rc').tolist()) - if vars_to_load: - vars_to_load = ComponentSet(vars_to_load) - iterator = filter(lambda var_rc: var_rc[0] in vars_to_load, iterator) - return ComponentMap(iterator) + res = ComponentMap() + gurobi_vars_to_load = [var_map[id(v)] for v in vars_to_load] + vals = solver_model.getAttr("Rc", gurobi_vars_to_load) + for var, val in zip(vars_to_load, vals): + res[var] = val -class GurobiSolverMixin: + return res + + +def _get_duals(solver_model, con_map, linear_cons_to_load, quadratic_cons_to_load): """ - gurobi_direct and gurobi_persistent check availability and set versions - in the same way. This moves the logic to a central location to reduce - duplicate code. + solver_model: gurobipy.Model + con_map: Dict[ConstraintData, gurobipy.Constr] + Maps the pyomo constraint to the gurobipy constraint + linear_cons_to_load: List[ConstraintData] + quadratic_cons_to_load: List[ConstraintData] """ + if solver_model.Status != gurobipy.GRB.OPTIMAL: + raise NoDualsError() + + linear_gurobi_cons = [con_map[c] for c in linear_cons_to_load] + quadratic_gurobi_cons = [con_map[c] for c in quadratic_cons_to_load] + linear_vals = solver_model.getAttr("Pi", linear_gurobi_cons) + quadratic_vals = solver_model.getAttr("QCPi", quadratic_gurobi_cons) + + duals = {} + for c, val in zip(linear_cons_to_load, linear_vals): + duals[c] = val + for c, val in zip(quadratic_cons_to_load, quadratic_vals): + duals[c] = val + return duals + + +class GurobiDirectBase(SolverBase): _num_gurobipy_env_clients = 0 _gurobipy_env = None _available = None _gurobipy_available = gurobipy_available + _tc_map = None + _minimum_version = (0, 0, 0) + + CONFIG = GurobiConfig() + + def __init__(self, **kwds): + super().__init__(**kwds) + self._register_env_client() + self._callback = None + + def __del__(self): + if not python_is_shutting_down(): + self._release_env_client() def available(self): if self._available is None: @@ -188,45 +218,47 @@ def available(self): # class* and not on the instance, or on the base class. That # allows different derived interfaces to have different # availability (e.g., persistent has a minimum version - # requirement that the direct interface doesn't) + # requirement that the direct interface doesn't - is that true?) if not self._gurobipy_available: if self._available is None: self.__class__._available = Availability.NotFound else: self.__class__._available = self._check_license() + if self.version() < self._minimum_version: + self.__class__._available = Availability.BadVersion return self._available @staticmethod def release_license(): - if GurobiSolverMixin._gurobipy_env is None: + if GurobiDirectBase._gurobipy_env is None: return - if GurobiSolverMixin._num_gurobipy_env_clients: + if GurobiDirectBase._num_gurobipy_env_clients: logger.warning( - "Call to GurobiSolverMixin.release_license() with %s remaining " - "environment clients." % (GurobiSolverMixin._num_gurobipy_env_clients,) + "Call to GurobiDirectBase.release_license() with %s remaining " + "environment clients." % (GurobiDirectBase._num_gurobipy_env_clients,) ) - GurobiSolverMixin._gurobipy_env.close() - GurobiSolverMixin._gurobipy_env = None + GurobiDirectBase._gurobipy_env.close() + GurobiDirectBase._gurobipy_env = None @staticmethod def env(): - if GurobiSolverMixin._gurobipy_env is None: + if GurobiDirectBase._gurobipy_env is None: with capture_output(capture_fd=True): - GurobiSolverMixin._gurobipy_env = gurobipy.Env() - return GurobiSolverMixin._gurobipy_env + GurobiDirectBase._gurobipy_env = gurobipy.Env() + return GurobiDirectBase._gurobipy_env @staticmethod def _register_env_client(): - GurobiSolverMixin._num_gurobipy_env_clients += 1 + GurobiDirectBase._num_gurobipy_env_clients += 1 @staticmethod def _release_env_client(): - GurobiSolverMixin._num_gurobipy_env_clients -= 1 - if GurobiSolverMixin._num_gurobipy_env_clients <= 0: + GurobiDirectBase._num_gurobipy_env_clients -= 1 + if GurobiDirectBase._num_gurobipy_env_clients <= 0: # Note that _num_gurobipy_env_clients should never be <0, # but if it is, release_license will issue a warning (that # we want to know about) - GurobiSolverMixin.release_license() + GurobiDirectBase.release_license() def _check_license(self): try: @@ -252,99 +284,49 @@ def version(self): ) return version + def _create_solver_model(self, pyomo_model): + # should return gurobi_model, solution_loader, has_objective + raise NotImplementedError('should be implemented by derived classes') -class GurobiDirect(GurobiSolverMixin, SolverBase): - """ - Interface to Gurobi using gurobipy - """ - - CONFIG = GurobiConfig() - - _tc_map = None - - def __init__(self, **kwds): - super().__init__(**kwds) - self._register_env_client() + def _pyomo_gurobi_var_iter(self): + # generator of tuples (pyomo_var, gurobi_var) + raise NotImplementedError('should be implemented by derived classes') - def __del__(self): - if not python_is_shutting_down(): - self._release_env_client() + def _mipstart(self): + for pyomo_var, gurobi_var in self._pyomo_gurobi_var_iter(): + if pyomo_var.is_integer() and pyomo_var.value is not None: + gurobi_var.setAttr('Start', pyomo_var.value) def solve(self, model, **kwds) -> Results: start_timestamp = datetime.datetime.now(datetime.timezone.utc) - config = self.config(value=kwds, preserve_implicit=True) - if not self.available(): - c = self.__class__ - raise ApplicationError( - f'Solver {c.__module__}.{c.__qualname__} is not available ' - f'({self.available()}).' - ) - if config.timer is None: - config.timer = HierarchicalTimer() - timer = config.timer - - StaleFlagManager.mark_all_as_stale() - - timer.start('compile_model') - repn = LinearStandardFormCompiler().write( - model, mixed_form=True, set_sense=None - ) - timer.stop('compile_model') + orig_config = self.config + orig_cwd = os.getcwd() + try: + config = self.config(value=kwds, preserve_implicit=True) - if len(repn.objectives) > 1: - raise IncompatibleModelError( - f"The {self.__class__.__name__} solver only supports models " - f"with zero or one objectives (received {len(repn.objectives)})." - ) + # hack to work around legacy solver wrapper __setattr__ + # otherwise, this would just be self.config = config + object.__setattr__(self, 'config', config) - timer.start('prepare_matrices') - inf = float('inf') - ninf = -inf - bounds = list(map(operator.attrgetter('bounds'), repn.columns)) - lb = [ninf if _b is None else _b for _b in map(operator.itemgetter(0), bounds)] - ub = [inf if _b is None else _b for _b in map(operator.itemgetter(1), bounds)] - CON = gurobipy.GRB.CONTINUOUS - BIN = gurobipy.GRB.BINARY - INT = gurobipy.GRB.INTEGER - vtype = [ - ( - CON - if v.is_continuous() - else BIN if v.is_binary() else INT if v.is_integer() else '?' - ) - for v in repn.columns - ] - sense_type = list('=<>') # Note: ordering matches 0, 1, -1 - sense = [sense_type[r[1]] for r in repn.rows] - timer.stop('prepare_matrices') + if not self.available(): + c = self.__class__ + raise ApplicationError( + f'Solver {c.__module__}.{c.__qualname__} is not available ' + f'({self.available()}).' + ) + if config.timer is None: + config.timer = HierarchicalTimer() + timer = config.timer - ostreams = [io.StringIO()] + config.tee - res = Results() + StaleFlagManager.mark_all_as_stale() + ostreams = [io.StringIO()] + config.tee - orig_cwd = os.getcwd() - try: if config.working_dir: os.chdir(config.working_dir) with capture_output(TeeStream(*ostreams), capture_fd=False): - gurobi_model = gurobipy.Model(env=self.env()) - - timer.start('transfer_model') - x = gurobi_model.addMVar( - len(repn.columns), - lb=lb, - ub=ub, - obj=repn.c.todense()[0] if repn.c.shape[0] else 0, - vtype=vtype, + gurobi_model, solution_loader, has_obj = self._create_solver_model( + model ) - A = gurobi_model.addMConstr(repn.A, x, sense, repn.rhs) - if repn.c.shape[0]: - gurobi_model.setAttr('ObjCon', repn.c_offset[0]) - gurobi_model.setAttr('ModelSense', int(repn.objectives[0].sense)) - # Note: calling gurobi_model.update() here is not - # necessary (it will happen as part of optimize()): - # gurobi_model.update() - timer.stop('transfer_model') - options = config.solver_options gurobi_model.setParam('LogToConsole', 1) @@ -359,42 +341,60 @@ def solve(self, model, **kwds) -> Results: gurobi_model.setParam('MIPGapAbs', config.abs_gap) if config.use_mipstart: - raise MouseTrap("MIPSTART not yet supported") + self._mipstart() for key, option in options.items(): gurobi_model.setParam(key, option) timer.start('optimize') - gurobi_model.optimize() + gurobi_model.optimize(self._callback) timer.stop('optimize') + + res = self._postsolve( + grb_model=gurobi_model, solution_loader=solution_loader, has_obj=has_obj + ) finally: os.chdir(orig_cwd) - res = self._postsolve( - timer, - config, - GurobiDirectSolutionLoader( - gurobi_model, A, x, repn.rows, repn.columns, repn.objectives - ), - ) + # hack to work around legacy solver wrapper __setattr__ + # otherwise, this would just be self.config = orig_config + object.__setattr__(self, 'config', orig_config) + self.config = orig_config - res.solver_config = config - res.solver_name = 'Gurobi' - res.solver_version = self.version() res.solver_log = ostreams[0].getvalue() - end_timestamp = datetime.datetime.now(datetime.timezone.utc) res.timing_info.start_timestamp = start_timestamp res.timing_info.wall_time = (end_timestamp - start_timestamp).total_seconds() res.timing_info.timer = timer return res - def _postsolve(self, timer: HierarchicalTimer, config, loader): - grb_model = loader._grb_model + def _get_tc_map(self): + if GurobiDirectBase._tc_map is None: + grb = gurobipy.GRB + tc = TerminationCondition + GurobiDirectBase._tc_map = { + grb.LOADED: tc.unknown, # problem is loaded, but no solution + grb.OPTIMAL: tc.convergenceCriteriaSatisfied, + grb.INFEASIBLE: tc.provenInfeasible, + grb.INF_OR_UNBD: tc.infeasibleOrUnbounded, + grb.UNBOUNDED: tc.unbounded, + grb.CUTOFF: tc.objectiveLimit, + grb.ITERATION_LIMIT: tc.iterationLimit, + grb.NODE_LIMIT: tc.iterationLimit, + grb.TIME_LIMIT: tc.maxTimeLimit, + grb.SOLUTION_LIMIT: tc.unknown, + grb.INTERRUPTED: tc.interrupted, + grb.NUMERIC: tc.unknown, + grb.SUBOPTIMAL: tc.unknown, + grb.USER_OBJ_LIMIT: tc.objectiveLimit, + } + return GurobiDirectBase._tc_map + + def _postsolve(self, grb_model, solution_loader, has_obj): status = grb_model.Status results = Results() - results.solution_loader = loader + results.solution_loader = solution_loader results.timing_info.gurobi_time = grb_model.Runtime if grb_model.SolCount > 0: @@ -412,11 +412,11 @@ def _postsolve(self, timer: HierarchicalTimer, config, loader): if ( results.termination_condition != TerminationCondition.convergenceCriteriaSatisfied - and config.raise_exception_on_nonoptimal_result + and self.config.raise_exception_on_nonoptimal_result ): raise NoOptimalSolutionError() - if loader._pyo_obj: + if has_obj: try: if math.isfinite(grb_model.ObjVal): results.incumbent_objective = grb_model.ObjVal @@ -437,34 +437,20 @@ def _postsolve(self, timer: HierarchicalTimer, config, loader): results.iteration_count = grb_model.getAttr('IterCount') - timer.start('load solution') - if config.load_solutions: + self.config.timer.start('load solution') + if self.config.load_solutions: if grb_model.SolCount > 0: results.solution_loader.load_vars() else: raise NoFeasibleSolutionError() - timer.stop('load solution') + self.config.timer.stop('load solution') - return results + # self.config gets copied a the beginning of + # solve and restored at the end, so modifying + # results.solver_config will not actually + # modify self.config + results.solver_config = self.config + results.solver_name = self.name + results.solver_version = self.version() - def _get_tc_map(self): - if GurobiDirect._tc_map is None: - grb = gurobipy.GRB - tc = TerminationCondition - GurobiDirect._tc_map = { - grb.LOADED: tc.unknown, # problem is loaded, but no solution - grb.OPTIMAL: tc.convergenceCriteriaSatisfied, - grb.INFEASIBLE: tc.provenInfeasible, - grb.INF_OR_UNBD: tc.infeasibleOrUnbounded, - grb.UNBOUNDED: tc.unbounded, - grb.CUTOFF: tc.objectiveLimit, - grb.ITERATION_LIMIT: tc.iterationLimit, - grb.NODE_LIMIT: tc.iterationLimit, - grb.TIME_LIMIT: tc.maxTimeLimit, - grb.SOLUTION_LIMIT: tc.unknown, - grb.INTERRUPTED: tc.interrupted, - grb.NUMERIC: tc.unknown, - grb.SUBOPTIMAL: tc.unknown, - grb.USER_OBJ_LIMIT: tc.objectiveLimit, - } - return GurobiDirect._tc_map + return results diff --git a/pyomo/contrib/solver/solvers/gurobi/gurobi_persistent.py b/pyomo/contrib/solver/solvers/gurobi/gurobi_persistent.py new file mode 100644 index 00000000000..603e8e21800 --- /dev/null +++ b/pyomo/contrib/solver/solvers/gurobi/gurobi_persistent.py @@ -0,0 +1,1385 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# 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 __future__ import annotations +import logging +from typing import Dict, List, Optional, Sequence, Mapping +from collections.abc import Iterable + +from pyomo.common.collections import ComponentSet, OrderedSet +from pyomo.common.shutdown import python_is_shutting_down +from pyomo.common.timing import HierarchicalTimer +from pyomo.core.base.objective import ObjectiveData +from pyomo.core.kernel.objective import minimize, maximize +from pyomo.core.base.var import VarData +from pyomo.core.base.constraint import ConstraintData, Constraint +from pyomo.core.base.sos import SOSConstraintData, SOSConstraint +from pyomo.core.base.param import ParamData +from pyomo.core.expr.numvalue import value, is_constant, is_fixed, native_numeric_types +from pyomo.repn import generate_standard_repn +from pyomo.contrib.solver.common.results import Results +from pyomo.contrib.solver.common.util import IncompatibleModelError +from pyomo.contrib.solver.common.solution_loader import SolutionLoaderBase +from pyomo.contrib.solver.common.base import PersistentSolverBase +from pyomo.core.staleflag import StaleFlagManager +from .gurobi_direct_base import ( + GurobiDirectBase, + gurobipy, + _load_vars, + _get_primals, + _get_duals, + _get_reduced_costs, +) +from pyomo.contrib.solver.common.util import get_objective +from pyomo.contrib.observer.model_observer import Observer, ModelChangeDetector + + +logger = logging.getLogger(__name__) + + +class GurobiDirectQuadraticSolutionLoader(SolutionLoaderBase): + def __init__( + self, solver_model, var_id_map, var_map, con_map, linear_cons, quadratic_cons + ) -> None: + super().__init__() + self._solver_model = solver_model + self._vars = var_id_map + self._var_map = var_map + self._con_map = con_map + self._linear_cons = linear_cons + self._quadratic_cons = quadratic_cons + GurobiDirectBase._register_env_client() + + def __del__(self): + # Release the gurobi license if this is the last reference to + # the environment (either through a results object or solver + # interface) + GurobiDirectBase._release_env_client() + + def load_vars( + self, vars_to_load: Optional[Sequence[VarData]] = None, solution_id=0 + ) -> None: + if vars_to_load is None: + vars_to_load = list(self._vars.values()) + _load_vars( + solver_model=self._solver_model, + var_map=self._var_map, + vars_to_load=vars_to_load, + solution_number=solution_id, + ) + + def get_primals( + self, vars_to_load: Optional[Sequence[VarData]] = None, solution_id=0 + ) -> Mapping[VarData, float]: + if vars_to_load is None: + vars_to_load = list(self._vars.values()) + return _get_primals( + solver_model=self._solver_model, + var_map=self._var_map, + vars_to_load=vars_to_load, + solution_number=solution_id, + ) + + def get_reduced_costs( + self, vars_to_load: Optional[Sequence[VarData]] = None + ) -> Mapping[VarData, float]: + if vars_to_load is None: + vars_to_load = list(self._vars.values()) + return _get_reduced_costs( + solver_model=self._solver_model, + var_map=self._var_map, + vars_to_load=vars_to_load, + ) + + def get_duals( + self, cons_to_load: Optional[Sequence[ConstraintData]] = None + ) -> Dict[ConstraintData, float]: + if cons_to_load is None: + cons_to_load = list(self._con_map.keys()) + linear_cons_to_load = [] + quadratic_cons_to_load = [] + for c in cons_to_load: + if c in self._linear_cons: + linear_cons_to_load.append(c) + else: + assert c in self._quadratic_cons + quadratic_cons_to_load.append(c) + return _get_duals( + solver_model=self._solver_model, + con_map=self._con_map, + linear_cons_to_load=linear_cons_to_load, + quadratic_cons_to_load=quadratic_cons_to_load, + ) + + +class GurobiPersistentSolutionLoader(GurobiDirectQuadraticSolutionLoader): + def __init__( + self, solver_model, var_id_map, var_map, con_map, linear_cons, quadratic_cons + ) -> None: + super().__init__( + solver_model, var_id_map, var_map, con_map, linear_cons, quadratic_cons + ) + self._valid = True + + def invalidate(self): + self._valid = False + + def _assert_solution_still_valid(self): + if not self._valid: + raise RuntimeError('The results in the solver are no longer valid.') + + def load_vars( + self, vars_to_load: Sequence[VarData] | None = None, solution_id=0 + ) -> None: + self._assert_solution_still_valid() + return super().load_vars(vars_to_load, solution_id) + + def get_primals( + self, vars_to_load: Sequence[VarData] | None = None, solution_id=0 + ) -> Mapping[VarData, float]: + self._assert_solution_still_valid() + return super().get_primals(vars_to_load, solution_id) + + def get_duals( + self, cons_to_load: Sequence[ConstraintData] | None = None + ) -> Dict[ConstraintData, float]: + self._assert_solution_still_valid() + return super().get_duals(cons_to_load) + + def get_reduced_costs( + self, vars_to_load: Sequence[VarData] | None = None + ) -> Mapping[VarData, float]: + self._assert_solution_still_valid() + return super().get_reduced_costs(vars_to_load) + + +class _MutableLowerBound: + def __init__(self, var_id, expr, var_map): + self.var_id = var_id + self.expr = expr + self.var_map = var_map + + def update(self): + self.var_map[self.var_id].setAttr('lb', value(self.expr)) + + +class _MutableUpperBound: + def __init__(self, var_id, expr, var_map): + self.var_id = var_id + self.expr = expr + self.var_map = var_map + + def update(self): + self.var_map[self.var_id].setAttr('ub', value(self.expr)) + + +class _MutableLinearCoefficient: + def __init__(self, expr, pyomo_con, con_map, pyomo_var_id, var_map, gurobi_model): + self.expr = expr + self.pyomo_con = pyomo_con + self.pyomo_var_id = pyomo_var_id + self.con_map = con_map + self.var_map = var_map + self.gurobi_model = gurobi_model + + @property + def gurobi_var(self): + return self.var_map[self.pyomo_var_id] + + @property + def gurobi_con(self): + return self.con_map[self.pyomo_con] + + def update(self): + self.gurobi_model.chgCoeff(self.gurobi_con, self.gurobi_var, value(self.expr)) + + +class _MutableRangeConstant: + def __init__( + self, lhs_expr, rhs_expr, pyomo_con, con_map, slack_name, gurobi_model + ): + self.lhs_expr = lhs_expr + self.rhs_expr = rhs_expr + self.pyomo_con = pyomo_con + self.con_map = con_map + self.slack_name = slack_name + self.gurobi_model = gurobi_model + + def update(self): + rhs_val = value(self.rhs_expr) + lhs_val = value(self.lhs_expr) + con = self.con_map[self.pyomo_con] + con.rhs = rhs_val + slack = self.gurobi_model.getVarByName(self.slack_name) + slack.ub = rhs_val - lhs_val + + +class _MutableConstant: + def __init__(self, expr, pyomo_con, con_map): + self.expr = expr + self.pyomo_con = pyomo_con + self.con_map = con_map + + def update(self): + con = self.con_map[self.pyomo_con] + con.rhs = value(self.expr) + + +class _MutableQuadraticConstraint: + def __init__( + self, gurobi_model, pyomo_con, con_map, constant, linear_coefs, quadratic_coefs + ): + self.pyomo_con = pyomo_con + self.con_map = con_map + self.gurobi_model = gurobi_model + self.constant = constant + self.last_constant_value = value(self.constant.expr) + self.linear_coefs = linear_coefs + self.last_linear_coef_values = [value(i.expr) for i in self.linear_coefs] + self.quadratic_coefs = quadratic_coefs + self.last_quadratic_coef_values = [value(i.expr) for i in self.quadratic_coefs] + + @property + def gurobi_con(self): + return self.con_map[self.pyomo_con] + + def get_updated_expression(self): + gurobi_expr = self.gurobi_model.getQCRow(self.gurobi_con) + for ndx, coef in enumerate(self.linear_coefs): + current_coef_value = value(coef.expr) + incremental_coef_value = ( + current_coef_value - self.last_linear_coef_values[ndx] + ) + gurobi_expr += incremental_coef_value * coef.gurobi_var + self.last_linear_coef_values[ndx] = current_coef_value + for ndx, coef in enumerate(self.quadratic_coefs): + current_coef_value = value(coef.expr) + incremental_coef_value = ( + current_coef_value - self.last_quadratic_coef_values[ndx] + ) + gurobi_expr += incremental_coef_value * coef.var1 * coef.var2 + self.last_quadratic_coef_values[ndx] = current_coef_value + return gurobi_expr + + def get_updated_rhs(self): + return value(self.constant.expr) + + +class _MutableObjective: + def __init__(self, gurobi_model, constant, linear_coefs, quadratic_coefs): + self.gurobi_model = gurobi_model + self.constant: _MutableConstant = constant + self.linear_coefs: List[_MutableLinearCoefficient] = linear_coefs + self.quadratic_coefs: List[_MutableQuadraticCoefficient] = quadratic_coefs + self.last_quadratic_coef_values: List[float] = [ + value(i.expr) for i in self.quadratic_coefs + ] + + def get_updated_expression(self): + for ndx, coef in enumerate(self.linear_coefs): + coef.gurobi_var.obj = value(coef.expr) + self.gurobi_model.ObjCon = value(self.constant.expr) + + gurobi_expr = None + for ndx, coef in enumerate(self.quadratic_coefs): + if value(coef.expr) != self.last_quadratic_coef_values[ndx]: + if gurobi_expr is None: + self.gurobi_model.update() + gurobi_expr = self.gurobi_model.getObjective() + current_coef_value = value(coef.expr) + incremental_coef_value = ( + current_coef_value - self.last_quadratic_coef_values[ndx] + ) + gurobi_expr += incremental_coef_value * coef.var1 * coef.var2 + self.last_quadratic_coef_values[ndx] = current_coef_value + return gurobi_expr + + +class _MutableQuadraticCoefficient: + def __init__(self, expr, v1id, v2id, var_map): + self.expr = expr + self.var_map = var_map + self.v1id = v1id + self.v2id = v2id + + @property + def var1(self): + return self.var_map[self.v1id] + + @property + def var2(self): + return self.var_map[self.v2id] + + +class GurobiDirectQuadratic(GurobiDirectBase): + _minimum_version = (7, 0, 0) + + def __init__(self, **kwds): + super().__init__(**kwds) + self._solver_model = None + self._vars = {} # from id(v) to v + self._pyomo_var_to_solver_var_map = {} + self._pyomo_con_to_solver_con_map = {} + self._linear_cons = set() + self._quadratic_cons = set() + self._pyomo_sos_to_solver_sos_map = {} + + def _create_solver_model(self, pyomo_model): + timer = self.config.timer + timer.start('create gurobipy model') + self._clear() + self._solver_model = gurobipy.Model(env=self.env()) + timer.start('collect constraints') + cons = list( + pyomo_model.component_data_objects( + Constraint, descend_into=True, active=True + ) + ) + timer.stop('collect constraints') + timer.start('translate constraints') + self._add_constraints(cons) + timer.stop('translate constraints') + timer.start('sos') + sos = list( + pyomo_model.component_data_objects( + SOSConstraint, descend_into=True, active=True + ) + ) + self._add_sos_constraints(sos) + timer.stop('sos') + timer.start('get objective') + obj = get_objective(pyomo_model) + timer.stop('get objective') + timer.start('translate objective') + self._set_objective(obj) + timer.stop('translate objective') + has_obj = obj is not None + solution_loader = GurobiDirectQuadraticSolutionLoader( + solver_model=self._solver_model, + var_id_map=self._vars, + var_map=self._pyomo_var_to_solver_var_map, + con_map=self._pyomo_con_to_solver_con_map, + linear_cons=self._linear_cons, + quadratic_cons=self._quadratic_cons, + ) + timer.stop('create gurobipy model') + return self._solver_model, solution_loader, has_obj + + def _clear(self): + self._solver_model = None + self._vars = {} + self._pyomo_var_to_solver_var_map = {} + self._pyomo_con_to_solver_con_map = {} + self._linear_cons = set() + self._quadratic_cons = set() + self._pyomo_sos_to_solver_sos_map = {} + + def _pyomo_gurobi_var_iter(self): + for vid, v in self._vars.items(): + yield v, self._pyomo_var_to_solver_var_map[vid] + + def _process_domain_and_bounds(self, var): + lb, ub, step = var.domain.get_interval() + if lb is None: + lb = -gurobipy.GRB.INFINITY + if ub is None: + ub = gurobipy.GRB.INFINITY + if step == 0: + vtype = gurobipy.GRB.CONTINUOUS + elif step == 1: + if lb == 0 and ub == 1: + vtype = gurobipy.GRB.BINARY + else: + vtype = gurobipy.GRB.INTEGER + else: + raise ValueError(f'Unrecognized domain: {var.domain}') + if var.fixed: + lb = var.value + ub = lb + else: + if var._lb is not None: + lb = max(lb, value(var._lb)) + if var._ub is not None: + ub = min(ub, value(var._ub)) + return lb, ub, vtype + + def _add_variables(self, variables: List[VarData]): + vtypes = [] + lbs = [] + ubs = [] + for ndx, var in enumerate(variables): + self._vars[id(var)] = var + lb, ub, vtype = self._process_domain_and_bounds(var) + vtypes.append(vtype) + lbs.append(lb) + ubs.append(ub) + + gurobi_vars = self._solver_model.addVars( + len(variables), lb=lbs, ub=ubs, vtype=vtypes + ).values() + + for pyomo_var, gurobi_var in zip(variables, gurobi_vars): + self._pyomo_var_to_solver_var_map[id(pyomo_var)] = gurobi_var + + def _get_expr_from_pyomo_repn(self, repn): + if repn.nonlinear_expr is not None: + raise IncompatibleModelError( + f'GurobiDirectQuadratic only supports linear and quadratic expressions: {repn}.' + ) + + if len(repn.linear_vars) > 0: + missing_vars = [v for v in repn.linear_vars if id(v) not in self._vars] + self._add_variables(missing_vars) + coef_list = [value(i) for i in repn.linear_coefs] + vlist = [self._pyomo_var_to_solver_var_map[id(v)] for v in repn.linear_vars] + new_expr = gurobipy.LinExpr(coef_list, vlist) + else: + new_expr = 0.0 + + if len(repn.quadratic_vars) > 0: + missing_vars = {} + for x, y in repn.quadratic_vars: + for v in [x, y]: + vid = id(v) + if vid not in self._vars: + missing_vars[vid] = v + self._add_variables(list(missing_vars.values())) + for coef, (x, y) in zip(repn.quadratic_coefs, repn.quadratic_vars): + gurobi_x = self._pyomo_var_to_solver_var_map[id(x)] + gurobi_y = self._pyomo_var_to_solver_var_map[id(y)] + new_expr += value(coef) * gurobi_x * gurobi_y + + return new_expr + + def _add_constraints(self, cons: List[ConstraintData]): + gurobi_expr_list = [] + for con in cons: + lb, body, ub = con.to_bounded_expression(evaluate_bounds=True) + repn = generate_standard_repn(body, quadratic=True, compute_values=True) + if len(repn.quadratic_vars) > 0: + self._quadratic_cons.add(con) + else: + self._linear_cons.add(con) + gurobi_expr = self._get_expr_from_pyomo_repn(repn) + if lb is None and ub is None: + raise ValueError( + "Constraint does not have a lower " f"or an upper bound: {con} \n" + ) + elif lb is None: + gurobi_expr_list.append(gurobi_expr <= float(ub - repn.constant)) + elif ub is None: + gurobi_expr_list.append(float(lb - repn.constant) <= gurobi_expr) + elif lb == ub: + gurobi_expr_list.append(gurobi_expr == float(lb - repn.constant)) + else: + gurobi_expr_list.append( + gurobi_expr + == [float(lb - repn.constant), float(ub - repn.constant)] + ) + + gurobi_cons = self._solver_model.addConstrs( + (gurobi_expr_list[i] for i in range(len(gurobi_expr_list))) + ).values() + self._pyomo_con_to_solver_con_map.update(zip(cons, gurobi_cons)) + + def _add_sos_constraints(self, cons: List[SOSConstraintData]): + for con in cons: + level = con.level + if level == 1: + sos_type = gurobipy.GRB.SOS_TYPE1 + elif level == 2: + sos_type = gurobipy.GRB.SOS_TYPE2 + else: + raise ValueError( + f"Solver does not support SOS level {level} constraints" + ) + + gurobi_vars = [] + weights = [] + + missing_vars = { + id(v): v for v, w in con.get_items() if id(v) not in self._vars + } + self._add_variables(list(missing_vars.values())) + + for v, w in con.get_items(): + v_id = id(v) + gurobi_vars.append(self._pyomo_var_to_solver_var_map[v_id]) + weights.append(w) + + gurobipy_con = self._solver_model.addSOS(sos_type, gurobi_vars, weights) + self._pyomo_sos_to_solver_sos_map[con] = gurobipy_con + + def _set_objective(self, obj): + if obj is None: + sense = gurobipy.GRB.MINIMIZE + gurobi_expr = 0 + repn_constant = 0 + else: + if obj.sense == minimize: + sense = gurobipy.GRB.MINIMIZE + elif obj.sense == maximize: + sense = gurobipy.GRB.MAXIMIZE + else: + raise ValueError(f'Objective sense is not recognized: {obj.sense}') + + repn = generate_standard_repn(obj.expr, quadratic=True, compute_values=True) + gurobi_expr = self._get_expr_from_pyomo_repn(repn) + repn_constant = repn.constant + + self._solver_model.setObjective(gurobi_expr + repn_constant, sense=sense) + + +class _GurobiObserver(Observer): + def __init__(self, opt: GurobiPersistentQuadratic) -> None: + self.opt = opt + + def add_variables(self, variables: List[VarData]): + self.opt._add_variables(variables) + + def add_parameters(self, params: List[ParamData]): + pass + + def add_constraints(self, cons: List[ConstraintData]): + self.opt._add_constraints(cons) + + def add_sos_constraints(self, cons: List[SOSConstraintData]): + self.opt._add_sos_constraints(cons) + + def set_objective(self, obj: ObjectiveData | None): + self.opt._set_objective(obj) + + def remove_constraints(self, cons: List[ConstraintData]): + self.opt._remove_constraints(cons) + + def remove_sos_constraints(self, cons: List[SOSConstraintData]): + self.opt._remove_sos_constraints(cons) + + def remove_variables(self, variables: List[VarData]): + self.opt._remove_variables(variables) + + def remove_parameters(self, params: List[ParamData]): + pass + + def update_variables(self, variables: List[VarData]): + self.opt._update_variables(variables) + + def update_parameters(self, params: List[ParamData]): + self.opt._update_parameters(params) + + +class GurobiPersistent(GurobiDirectQuadratic, PersistentSolverBase): + _minimum_version = (7, 0, 0) + + def __init__(self, **kwds): + super().__init__(**kwds) + self._pyomo_model = None + self._objective = None + self._mutable_helpers = {} + self._mutable_bounds = {} + self._mutable_quadratic_helpers = {} + self._mutable_objective = None + self._needs_updated = True + self._callback_func = None + self._constraints_added_since_update = OrderedSet() + self._vars_added_since_update = ComponentSet() + self._last_results_object: Optional[Results] = None + self._observer = _GurobiObserver(self) + self._change_detector = ModelChangeDetector(observers=[self._observer]) + self._constraint_ndx = 0 + self._should_update_parameters = False + + @property + def auto_updates(self): + return self._change_detector.config + + def _clear(self): + super()._clear() + self._pyomo_model = None + self._objective = None + self._mutable_helpers = {} + self._mutable_bounds = {} + self._mutable_quadratic_helpers = {} + self._mutable_objective = None + self._needs_updated = True + self._constraints_added_since_update = OrderedSet() + self._vars_added_since_update = ComponentSet() + self._last_results_object = None + self._constraint_ndx = 0 + + def _create_solver_model(self, pyomo_model): + if pyomo_model is self._pyomo_model: + self.update() + else: + self.set_instance(pyomo_model) + + solution_loader = GurobiPersistentSolutionLoader( + solver_model=self._solver_model, + var_id_map=self._vars, + var_map=self._pyomo_var_to_solver_var_map, + con_map=self._pyomo_con_to_solver_con_map, + linear_cons=self._linear_cons, + quadratic_cons=self._quadratic_cons, + ) + has_obj = self._objective is not None + return self._solver_model, solution_loader, has_obj + + def release_license(self): + self._clear() + self.__class__.release_license() + + def solve(self, model, **kwds) -> Results: + res = super().solve(model, **kwds) + self._needs_updated = False + return res + + def _process_domain_and_bounds(self, var): + res = super()._process_domain_and_bounds(var) + if not is_constant(var._lb): + mutable_lb = _MutableLowerBound( + id(var), var.lower, self._pyomo_var_to_solver_var_map + ) + self._mutable_bounds[id(var), 'lb'] = (var, mutable_lb) + if not is_constant(var._ub): + mutable_ub = _MutableUpperBound( + id(var), var.upper, self._pyomo_var_to_solver_var_map + ) + self._mutable_bounds[id(var), 'ub'] = (var, mutable_ub) + return res + + def _add_variables(self, variables: List[VarData]): + self._invalidate_last_results() + super()._add_variables(variables) + self._vars_added_since_update.update(variables) + self._needs_updated = True + + def set_instance(self, pyomo_model): + if self.config.timer is None: + timer = HierarchicalTimer() + else: + timer = self.config.timer + self._clear() + self._pyomo_model = pyomo_model + self._solver_model = gurobipy.Model(env=self.env()) + timer.start('set_instance') + self._change_detector.set_instance(pyomo_model) + timer.stop('set_instance') + + def update(self): + if self.config.timer is None: + timer = HierarchicalTimer() + else: + timer = self.config.timer + if self._pyomo_model is None: + raise RuntimeError('must call set_instance or solve before update') + timer.start('update') + if self._needs_updated: + self._update_gurobi_model() + self._change_detector.update(timer=timer) + if self._should_update_parameters: + self._update_parameters([]) + timer.stop('update') + + def _add_constraints(self, cons: List[ConstraintData]): + self._invalidate_last_results() + gurobi_expr_list = [] + for ndx, con in enumerate(cons): + lb, body, ub = con.to_bounded_expression(evaluate_bounds=False) + repn = generate_standard_repn(body, quadratic=True, compute_values=False) + if len(repn.quadratic_vars) > 0: + self._quadratic_cons.add(con) + else: + self._linear_cons.add(con) + gurobi_expr = self._get_expr_from_pyomo_repn(repn) + mutable_constant = None + if lb is None and ub is None: + raise ValueError( + "Constraint does not have a lower " f"or an upper bound: {con} \n" + ) + elif lb is None: + rhs_expr = ub - repn.constant + gurobi_expr_list.append(gurobi_expr <= float(value(rhs_expr))) + if not is_constant(rhs_expr): + mutable_constant = _MutableConstant( + rhs_expr, con, self._pyomo_con_to_solver_con_map + ) + elif ub is None: + rhs_expr = lb - repn.constant + gurobi_expr_list.append(float(value(rhs_expr)) <= gurobi_expr) + if not is_constant(rhs_expr): + mutable_constant = _MutableConstant( + rhs_expr, con, self._pyomo_con_to_solver_con_map + ) + elif con.equality: + rhs_expr = lb - repn.constant + gurobi_expr_list.append(gurobi_expr == float(value(rhs_expr))) + if not is_constant(rhs_expr): + mutable_constant = _MutableConstant( + rhs_expr, con, self._pyomo_con_to_solver_con_map + ) + else: + assert ( + len(repn.quadratic_vars) == 0 + ), "Quadratic range constraints are not supported" + lhs_expr = lb - repn.constant + rhs_expr = ub - repn.constant + gurobi_expr_list.append( + gurobi_expr == [float(value(lhs_expr)), float(value(rhs_expr))] + ) + if not is_constant(lhs_expr) or not is_constant(rhs_expr): + conname = f'c{self._constraint_ndx}[{ndx}]' + mutable_constant = _MutableRangeConstant( + lhs_expr, + rhs_expr, + con, + self._pyomo_con_to_solver_con_map, + 'Rg' + conname, + self._solver_model, + ) + + mlc_list = [] + for c, v in zip(repn.linear_coefs, repn.linear_vars): + if not is_constant(c): + mlc = _MutableLinearCoefficient( + c, + con, + self._pyomo_con_to_solver_con_map, + id(v), + self._pyomo_var_to_solver_var_map, + self._solver_model, + ) + mlc_list.append(mlc) + + if len(repn.quadratic_vars) == 0: + if len(mlc_list) > 0: + self._mutable_helpers[con] = mlc_list + if mutable_constant is not None: + if con not in self._mutable_helpers: + self._mutable_helpers[con] = [] + self._mutable_helpers[con].append(mutable_constant) + else: + if mutable_constant is None: + mutable_constant = _MutableConstant( + rhs_expr, con, self._pyomo_con_to_solver_con_map + ) + mqc_list = [] + for coef, (x, y) in zip(repn.quadratic_coefs, repn.quadratic_vars): + if not is_constant(coef): + mqc = _MutableQuadraticCoefficient( + coef, id(x), id(y), self._pyomo_var_to_solver_var_map + ) + mqc_list.append(mqc) + mqc = _MutableQuadraticConstraint( + self._solver_model, + con, + self._pyomo_con_to_solver_con_map, + mutable_constant, + mlc_list, + mqc_list, + ) + self._mutable_quadratic_helpers[con] = mqc + + gurobi_cons = list( + self._solver_model.addConstrs( + (gurobi_expr_list[i] for i in range(len(gurobi_expr_list))), + name=f'c{self._constraint_ndx}', + ).values() + ) + self._constraint_ndx += 1 + self._pyomo_con_to_solver_con_map.update(zip(cons, gurobi_cons)) + self._constraints_added_since_update.update(cons) + self._needs_updated = True + + def _add_sos_constraints(self, cons: List[SOSConstraintData]): + self._invalidate_last_results() + super()._add_sos_constraints(cons) + self._constraints_added_since_update.update(cons) + self._needs_updated = True + + def _set_objective(self, obj): + self._invalidate_last_results() + if obj is None: + sense = gurobipy.GRB.MINIMIZE + gurobi_expr = 0 + repn_constant = 0 + self._mutable_objective = None + else: + if obj.sense == minimize: + sense = gurobipy.GRB.MINIMIZE + elif obj.sense == maximize: + sense = gurobipy.GRB.MAXIMIZE + else: + raise ValueError(f'Objective sense is not recognized: {obj.sense}') + + repn = generate_standard_repn( + obj.expr, quadratic=True, compute_values=False + ) + repn_constant = value(repn.constant) + gurobi_expr = self._get_expr_from_pyomo_repn(repn) + + mutable_constant = _MutableConstant(repn.constant, None, None) + + mlc_list = [] + for c, v in zip(repn.linear_coefs, repn.linear_vars): + if not is_constant(c): + mlc = _MutableLinearCoefficient( + c, + None, + None, + id(v), + self._pyomo_var_to_solver_var_map, + self._solver_model, + ) + mlc_list.append(mlc) + + mqc_list = [] + for coef, (x, y) in zip(repn.quadratic_coefs, repn.quadratic_vars): + if not is_constant(coef): + mqc = _MutableQuadraticCoefficient( + coef, id(x), id(y), self._pyomo_var_to_solver_var_map + ) + mqc_list.append(mqc) + + self._mutable_objective = _MutableObjective( + self._solver_model, mutable_constant, mlc_list, mqc_list + ) + + # hack + # see PR #2454 + if self._objective is not None: + self._solver_model.setObjective(0) + self._solver_model.update() + + self._solver_model.setObjective(gurobi_expr + repn_constant, sense=sense) + self._objective = obj + self._needs_updated = True + + def _update_gurobi_model(self): + self._solver_model.update() + self._constraints_added_since_update = OrderedSet() + self._vars_added_since_update = ComponentSet() + self._needs_updated = False + + def _remove_constraints(self, cons: List[ConstraintData]): + self._invalidate_last_results() + for con in cons: + if con in self._constraints_added_since_update: + self._update_gurobi_model() + solver_con = self._pyomo_con_to_solver_con_map[con] + self._solver_model.remove(solver_con) + del self._pyomo_con_to_solver_con_map[con] + self._mutable_helpers.pop(con, None) + self._mutable_quadratic_helpers.pop(con, None) + self._needs_updated = True + + def _remove_sos_constraints(self, cons: List[SOSConstraintData]): + self._invalidate_last_results() + for con in cons: + if con in self._constraints_added_since_update: + self._update_gurobi_model() + solver_sos_con = self._pyomo_sos_to_solver_sos_map[con] + self._solver_model.remove(solver_sos_con) + del self._pyomo_sos_to_solver_sos_map[con] + self._needs_updated = True + + def _remove_variables(self, variables: List[VarData]): + self._invalidate_last_results() + for var in variables: + v_id = id(var) + if var in self._vars_added_since_update: + self._update_gurobi_model() + solver_var = self._pyomo_var_to_solver_var_map[v_id] + self._solver_model.remove(solver_var) + del self._pyomo_var_to_solver_var_map[v_id] + del self._vars[v_id] + self._mutable_bounds.pop(v_id, None) + self._needs_updated = True + + def _update_variables(self, variables: List[VarData]): + self._invalidate_last_results() + for var in variables: + var_id = id(var) + if var_id not in self._pyomo_var_to_solver_var_map: + raise ValueError( + f'The Var provided to update_var needs to be added first: {var}' + ) + self._mutable_bounds.pop((var_id, 'lb'), None) + self._mutable_bounds.pop((var_id, 'ub'), None) + gurobipy_var = self._pyomo_var_to_solver_var_map[var_id] + lb, ub, vtype = self._process_domain_and_bounds(var) + gurobipy_var.setAttr('lb', lb) + gurobipy_var.setAttr('ub', ub) + gurobipy_var.setAttr('vtype', vtype) + if var.fixed: + self._should_update_parameters = True + self._needs_updated = True + + def _update_parameters(self, params: List[ParamData]): + self._invalidate_last_results() + for con, helpers in self._mutable_helpers.items(): + for helper in helpers: + helper.update() + for k, (v, helper) in self._mutable_bounds.items(): + helper.update() + + for con, helper in self._mutable_quadratic_helpers.items(): + if con in self._constraints_added_since_update: + self._update_gurobi_model() + gurobi_con = helper.gurobi_con + new_gurobi_expr = helper.get_updated_expression() + new_rhs = helper.get_updated_rhs() + new_sense = gurobi_con.qcsense + self._solver_model.remove(gurobi_con) + new_con = self._solver_model.addQConstr(new_gurobi_expr, new_sense, new_rhs) + self._pyomo_con_to_solver_con_map[con] = new_con + helper.pyomo_con = con + self._constraints_added_since_update.add(con) + + if self._mutable_objective is not None: + new_gurobi_expr = self._mutable_objective.get_updated_expression() + if new_gurobi_expr is not None: + if self._objective.sense == minimize: + sense = gurobipy.GRB.MINIMIZE + else: + sense = gurobipy.GRB.MAXIMIZE + # TODO: need a test for when part of the object is linear + # and part of the objective is quadratic, but both + # parts have mutable coefficients + self._solver_model.setObjective(new_gurobi_expr, sense=sense) + + self._should_update_parameters = False + + def _invalidate_last_results(self): + if self._last_results_object is not None: + self._last_results_object.solution_loader.invalidate() + + def get_model_attr(self, attr): + """ + Get the value of an attribute on the Gurobi model. + + Parameters + ---------- + attr: str + The attribute to get. See Gurobi documentation for descriptions of the attributes. + """ + if self._needs_updated: + self._update_gurobi_model() + return self._solver_model.getAttr(attr) + + def write(self, filename): + """ + Write the model to a file (e.g., and lp file). + + Parameters + ---------- + filename: str + Name of the file to which the model should be written. + """ + self._solver_model.write(filename) + self._constraints_added_since_update = OrderedSet() + self._vars_added_since_update = ComponentSet() + self._needs_updated = False + + def set_linear_constraint_attr(self, con, attr, val): + """ + Set the value of an attribute on a gurobi linear constraint. + + Parameters + ---------- + con: pyomo.core.base.constraint.ConstraintData + The pyomo constraint for which the corresponding gurobi constraint attribute + should be modified. + attr: str + The attribute to be modified. Options are: + CBasis + DStart + Lazy + val: any + See gurobi documentation for acceptable values. + """ + if attr in {'Sense', 'RHS', 'ConstrName'}: + raise ValueError( + f'Linear constraint attr {attr} cannot be set with' + ' the set_linear_constraint_attr method. Please use' + ' the remove_constraint and add_constraint methods.' + ) + self._pyomo_con_to_solver_con_map[con].setAttr(attr, val) + self._needs_updated = True + + def set_var_attr(self, var, attr, val): + """ + Set the value of an attribute on a gurobi variable. + + Parameters + ---------- + var: pyomo.core.base.var.VarData + The pyomo var for which the corresponding gurobi var attribute + should be modified. + attr: str + The attribute to be modified. Options are: + Start + VarHintVal + VarHintPri + BranchPriority + VBasis + PStart + val: any + See gurobi documentation for acceptable values. + """ + if attr in {'LB', 'UB', 'VType', 'VarName'}: + raise ValueError( + f'Var attr {attr} cannot be set with' + ' the set_var_attr method. Please use' + ' the update_var method.' + ) + if attr == 'Obj': + raise ValueError( + 'Var attr Obj cannot be set with' + ' the set_var_attr method. Please use' + ' the set_objective method.' + ) + self._pyomo_var_to_solver_var_map[id(var)].setAttr(attr, val) + self._needs_updated = True + + def get_var_attr(self, var, attr): + """ + Get the value of an attribute on a gurobi var. + + Parameters + ---------- + var: pyomo.core.base.var.VarData + The pyomo var for which the corresponding gurobi var attribute + should be retrieved. + attr: str + The attribute to get. See gurobi documentation + """ + if self._needs_updated: + self._update_gurobi_model() + return self._pyomo_var_to_solver_var_map[id(var)].getAttr(attr) + + def get_linear_constraint_attr(self, con, attr): + """ + Get the value of an attribute on a gurobi linear constraint. + + Parameters + ---------- + con: pyomo.core.base.constraint.ConstraintData + The pyomo constraint for which the corresponding gurobi constraint attribute + should be retrieved. + attr: str + The attribute to get. See the Gurobi documentation + """ + if self._needs_updated: + self._update_gurobi_model() + return self._pyomo_con_to_solver_con_map[con].getAttr(attr) + + def get_sos_attr(self, con, attr): + """ + Get the value of an attribute on a gurobi sos constraint. + + Parameters + ---------- + con: pyomo.core.base.sos.SOSConstraintData + The pyomo SOS constraint for which the corresponding gurobi SOS constraint attribute + should be retrieved. + attr: str + The attribute to get. See the Gurobi documentation + """ + if self._needs_updated: + self._update_gurobi_model() + return self._pyomo_sos_to_solver_sos_map[con].getAttr(attr) + + def get_quadratic_constraint_attr(self, con, attr): + """ + Get the value of an attribute on a gurobi quadratic constraint. + + Parameters + ---------- + con: pyomo.core.base.constraint.ConstraintData + The pyomo constraint for which the corresponding gurobi constraint attribute + should be retrieved. + attr: str + The attribute to get. See the Gurobi documentation + """ + if self._needs_updated: + self._update_gurobi_model() + return self._pyomo_con_to_solver_con_map[con].getAttr(attr) + + def set_gurobi_param(self, param, val): + """ + Set a gurobi parameter. + + Parameters + ---------- + param: str + The gurobi parameter to set. Options include any gurobi parameter. + Please see the Gurobi documentation for options. + val: any + The value to set the parameter to. See Gurobi documentation for possible values. + """ + self._solver_model.setParam(param, val) + + def get_gurobi_param_info(self, param): + """ + Get information about a gurobi parameter. + + Parameters + ---------- + param: str + The gurobi parameter to get info for. See Gurobi documentation for possible options. + + Returns + ------- + six-tuple containing the parameter name, type, value, minimum value, maximum value, and default value. + """ + return self._solver_model.getParamInfo(param) + + def _intermediate_callback(self): + def f(gurobi_model, where): + self._callback_func(self._pyomo_model, self, where) + + return f + + def set_callback(self, func=None): + """ + Specify a callback for gurobi to use. + + Parameters + ---------- + func: function + The function to call. The function should have three arguments. The first will be the pyomo model being + solved. The second will be the GurobiPersistent instance. The third will be an enum member of + gurobipy.GRB.Callback. This will indicate where in the branch and bound algorithm gurobi is at. For + example, suppose we want to solve + + .. math:: + + min 2*x + y + + s.t. + + y >= (x-2)**2 + + 0 <= x <= 4 + + y >= 0 + + y integer + + as an MILP using extended cutting planes in callbacks. + + >>> from gurobipy import GRB # doctest:+SKIP + >>> import pyomo.environ as pyo + >>> from pyomo.core.expr.taylor_series import taylor_series_expansion + >>> from pyomo.contrib import appsi + >>> + >>> m = pyo.ConcreteModel() + >>> m.x = pyo.Var(bounds=(0, 4)) + >>> m.y = pyo.Var(within=pyo.Integers, bounds=(0, None)) + >>> m.obj = pyo.Objective(expr=2*m.x + m.y) + >>> m.cons = pyo.ConstraintList() # for the cutting planes + >>> + >>> def _add_cut(xval): + ... # a function to generate the cut + ... m.x.value = xval + ... return m.cons.add(m.y >= taylor_series_expansion((m.x - 2)**2)) + ... + >>> _c = _add_cut(0) # start with 2 cuts at the bounds of x + >>> _c = _add_cut(4) # this is an arbitrary choice + >>> + >>> opt = appsi.solvers.Gurobi() + >>> opt.config.stream_solver = True + >>> opt.set_instance(m) # doctest:+SKIP + >>> opt.gurobi_options['PreCrush'] = 1 + >>> opt.gurobi_options['LazyConstraints'] = 1 + >>> + >>> def my_callback(cb_m, cb_opt, cb_where): + ... if cb_where == GRB.Callback.MIPSOL: + ... cb_opt.cbGetSolution(variables=[m.x, m.y]) + ... if m.y.value < (m.x.value - 2)**2 - 1e-6: + ... cb_opt.cbLazy(_add_cut(m.x.value)) + ... + >>> opt.set_callback(my_callback) + >>> res = opt.solve(m) # doctest:+SKIP + + """ + if func is not None: + self._callback_func = func + self._callback = self._intermediate_callback() + else: + self._callback = None + self._callback_func = None + + def cbCut(self, con): + """ + Add a cut within a callback. + + Parameters + ---------- + con: pyomo.core.base.constraint.ConstraintData + The cut to add + """ + if not con.active: + raise ValueError('cbCut expected an active constraint.') + + if is_fixed(con.body): + raise ValueError('cbCut expected a non-trivial constraint') + + repn = generate_standard_repn(con.body, quadratic=True, compute_values=True) + gurobi_expr = self._get_expr_from_pyomo_repn(repn) + + if con.has_lb(): + if con.has_ub(): + raise ValueError('Range constraints are not supported in cbCut.') + if not is_fixed(con.lower): + raise ValueError(f'Lower bound of constraint {con} is not constant.') + if con.has_ub(): + if not is_fixed(con.upper): + raise ValueError(f'Upper bound of constraint {con} is not constant.') + + if con.equality: + self._solver_model.cbCut( + lhs=gurobi_expr, + sense=gurobipy.GRB.EQUAL, + rhs=value(con.lower - repn.constant), + ) + elif con.has_lb() and (value(con.lower) > -float('inf')): + self._solver_model.cbCut( + lhs=gurobi_expr, + sense=gurobipy.GRB.GREATER_EQUAL, + rhs=value(con.lower - repn.constant), + ) + elif con.has_ub() and (value(con.upper) < float('inf')): + self._solver_model.cbCut( + lhs=gurobi_expr, + sense=gurobipy.GRB.LESS_EQUAL, + rhs=value(con.upper - repn.constant), + ) + else: + raise ValueError( + f'Constraint does not have a lower or an upper bound {con} \n' + ) + + def cbGet(self, what): + return self._solver_model.cbGet(what) + + def cbGetNodeRel(self, variables): + """ + Parameters + ---------- + variables: Var or iterable of Var + """ + if not isinstance(variables, Iterable): + variables = [variables] + gurobi_vars = [self._pyomo_var_to_solver_var_map[id(i)] for i in variables] + var_values = self._solver_model.cbGetNodeRel(gurobi_vars) + for i, v in enumerate(variables): + v.set_value(var_values[i], skip_validation=True) + + def cbGetSolution(self, variables): + """ + Parameters + ---------- + variables: iterable of vars + """ + if not isinstance(variables, Iterable): + variables = [variables] + gurobi_vars = [self._pyomo_var_to_solver_var_map[id(i)] for i in variables] + var_values = self._solver_model.cbGetSolution(gurobi_vars) + for i, v in enumerate(variables): + v.set_value(var_values[i], skip_validation=True) + + def cbLazy(self, con): + """ + Parameters + ---------- + con: pyomo.core.base.constraint.ConstraintData + The lazy constraint to add + """ + if not con.active: + raise ValueError('cbLazy expected an active constraint.') + + if is_fixed(con.body): + raise ValueError('cbLazy expected a non-trivial constraint') + + repn = generate_standard_repn(con.body, quadratic=True, compute_values=True) + gurobi_expr = self._get_expr_from_pyomo_repn(repn) + + if con.has_lb(): + if con.has_ub(): + raise ValueError('Range constraints are not supported in cbLazy.') + if not is_fixed(con.lower): + raise ValueError(f'Lower bound of constraint {con} is not constant.') + if con.has_ub(): + if not is_fixed(con.upper): + raise ValueError(f'Upper bound of constraint {con} is not constant.') + + if con.equality: + self._solver_model.cbLazy( + lhs=gurobi_expr, + sense=gurobipy.GRB.EQUAL, + rhs=value(con.lower - repn.constant), + ) + elif con.has_lb() and (value(con.lower) > -float('inf')): + self._solver_model.cbLazy( + lhs=gurobi_expr, + sense=gurobipy.GRB.GREATER_EQUAL, + rhs=value(con.lower - repn.constant), + ) + elif con.has_ub() and (value(con.upper) < float('inf')): + self._solver_model.cbLazy( + lhs=gurobi_expr, + sense=gurobipy.GRB.LESS_EQUAL, + rhs=value(con.upper - repn.constant), + ) + else: + raise ValueError( + f'Constraint does not have a lower or an upper bound {con} \n' + ) + + def cbSetSolution(self, variables, solution): + if not isinstance(variables, Iterable): + variables = [variables] + gurobi_vars = [self._pyomo_var_to_solver_var_map[id(i)] for i in variables] + self._solver_model.cbSetSolution(gurobi_vars, solution) + + def cbUseSolution(self): + return self._solver_model.cbUseSolution() + + def reset(self): + self._solver_model.reset() + + def add_variables(self, variables): + self._change_detector.add_variables(variables) + + def add_constraints(self, cons): + self._change_detector.add_constraints(cons) + + def add_sos_constraints(self, cons): + self._change_detector.add_sos_constraints(cons) + + def set_objective(self, obj): + self._change_detector.set_objective(obj) + + def remove_constraints(self, cons): + self._change_detector.remove_constraints(cons) + + def remove_sos_constraints(self, cons): + self._change_detector.remove_sos_constraints(cons) + + def remove_variables(self, variables): + self._change_detector.remove_variables(variables) + + def update_variables(self, variables): + self._change_detector.update_variables(variables) + + def update_parameters(self, params): + self._change_detector.update_parameters(params) diff --git a/pyomo/contrib/solver/solvers/gurobi_persistent.py b/pyomo/contrib/solver/solvers/gurobi_persistent.py deleted file mode 100644 index ea3693c1c70..00000000000 --- a/pyomo/contrib/solver/solvers/gurobi_persistent.py +++ /dev/null @@ -1,1409 +0,0 @@ -# ___________________________________________________________________________ -# -# Pyomo: Python Optimization Modeling Objects -# Copyright (c) 2008-2025 -# 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. -# ___________________________________________________________________________ - -import io -import logging -import math -from typing import List, Optional -from collections.abc import Iterable - -from pyomo.common.collections import ComponentSet, ComponentMap, OrderedSet -from pyomo.common.dependencies import attempt_import -from pyomo.common.errors import ApplicationError -from pyomo.common.tee import capture_output, TeeStream -from pyomo.common.timing import HierarchicalTimer -from pyomo.common.shutdown import python_is_shutting_down -from pyomo.core.kernel.objective import minimize, maximize -from pyomo.core.base import SymbolMap, NumericLabeler, TextLabeler -from pyomo.core.base.var import VarData -from pyomo.core.base.constraint import ConstraintData -from pyomo.core.base.sos import SOSConstraintData -from pyomo.core.base.param import ParamData -from pyomo.core.expr.numvalue import value, is_constant, is_fixed, native_numeric_types -from pyomo.repn import generate_standard_repn -from pyomo.core.expr.numeric_expr import NPV_MaxExpression, NPV_MinExpression -from pyomo.contrib.solver.common.base import PersistentSolverBase, Availability -from pyomo.contrib.solver.common.results import ( - Results, - TerminationCondition, - SolutionStatus, -) -from pyomo.contrib.solver.common.config import PersistentBranchAndBoundConfig -from pyomo.contrib.solver.solvers.gurobi_direct import ( - GurobiConfigMixin, - GurobiSolverMixin, -) -from pyomo.contrib.solver.common.util import ( - NoFeasibleSolutionError, - NoOptimalSolutionError, - NoDualsError, - NoReducedCostsError, - NoSolutionError, - IncompatibleModelError, -) -from pyomo.contrib.solver.common.persistent import ( - PersistentSolverUtils, - PersistentSolverMixin, -) -from pyomo.contrib.solver.common.solution_loader import PersistentSolutionLoader -from pyomo.core.staleflag import StaleFlagManager - - -logger = logging.getLogger(__name__) - - -def _import_gurobipy(): - try: - import gurobipy - except ImportError: - GurobiPersistent._available = Availability.NotFound - raise - if gurobipy.GRB.VERSION_MAJOR < 7: - GurobiPersistent._available = Availability.BadVersion - raise ImportError('The Persistent Gurobi interface requires gurobipy>=7.0.0') - return gurobipy - - -gurobipy, gurobipy_available = attempt_import('gurobipy', importer=_import_gurobipy) - - -class GurobiConfig(PersistentBranchAndBoundConfig, GurobiConfigMixin): - def __init__( - self, - description=None, - doc=None, - implicit=False, - implicit_domain=None, - visibility=0, - ): - PersistentBranchAndBoundConfig.__init__( - self, - description=description, - doc=doc, - implicit=implicit, - implicit_domain=implicit_domain, - visibility=visibility, - ) - GurobiConfigMixin.__init__(self) - - -class GurobiSolutionLoader(PersistentSolutionLoader): - def load_vars(self, vars_to_load=None, solution_number=0): - self._assert_solution_still_valid() - self._solver._load_vars( - vars_to_load=vars_to_load, solution_number=solution_number - ) - - def get_primals(self, vars_to_load=None, solution_number=0): - self._assert_solution_still_valid() - return self._solver._get_primals( - vars_to_load=vars_to_load, solution_number=solution_number - ) - - -class _MutableLowerBound: - def __init__(self, expr): - self.var = None - self.expr = expr - - def update(self): - self.var.setAttr('lb', value(self.expr)) - - -class _MutableUpperBound: - def __init__(self, expr): - self.var = None - self.expr = expr - - def update(self): - self.var.setAttr('ub', value(self.expr)) - - -class _MutableLinearCoefficient: - def __init__(self): - self.expr = None - self.var = None - self.con = None - self.gurobi_model = None - - def update(self): - self.gurobi_model.chgCoeff(self.con, self.var, value(self.expr)) - - -class _MutableRangeConstant: - def __init__(self): - self.lhs_expr = None - self.rhs_expr = None - self.con = None - self.slack_name = None - self.gurobi_model = None - - def update(self): - rhs_val = value(self.rhs_expr) - lhs_val = value(self.lhs_expr) - self.con.rhs = rhs_val - slack = self.gurobi_model.getVarByName(self.slack_name) - slack.ub = rhs_val - lhs_val - - -class _MutableConstant: - def __init__(self): - self.expr = None - self.con = None - - def update(self): - self.con.rhs = value(self.expr) - - -class _MutableQuadraticConstraint: - def __init__( - self, gurobi_model, gurobi_con, constant, linear_coefs, quadratic_coefs - ): - self.con = gurobi_con - self.gurobi_model = gurobi_model - self.constant = constant - self.last_constant_value = value(self.constant.expr) - self.linear_coefs = linear_coefs - self.last_linear_coef_values = [value(i.expr) for i in self.linear_coefs] - self.quadratic_coefs = quadratic_coefs - self.last_quadratic_coef_values = [value(i.expr) for i in self.quadratic_coefs] - - def get_updated_expression(self): - gurobi_expr = self.gurobi_model.getQCRow(self.con) - for ndx, coef in enumerate(self.linear_coefs): - current_coef_value = value(coef.expr) - incremental_coef_value = ( - current_coef_value - self.last_linear_coef_values[ndx] - ) - gurobi_expr += incremental_coef_value * coef.var - self.last_linear_coef_values[ndx] = current_coef_value - for ndx, coef in enumerate(self.quadratic_coefs): - current_coef_value = value(coef.expr) - incremental_coef_value = ( - current_coef_value - self.last_quadratic_coef_values[ndx] - ) - gurobi_expr += incremental_coef_value * coef.var1 * coef.var2 - self.last_quadratic_coef_values[ndx] = current_coef_value - return gurobi_expr - - def get_updated_rhs(self): - return value(self.constant.expr) - - -class _MutableObjective: - def __init__(self, gurobi_model, constant, linear_coefs, quadratic_coefs): - self.gurobi_model = gurobi_model - self.constant = constant - self.linear_coefs = linear_coefs - self.quadratic_coefs = quadratic_coefs - self.last_quadratic_coef_values = [value(i.expr) for i in self.quadratic_coefs] - - def get_updated_expression(self): - for ndx, coef in enumerate(self.linear_coefs): - coef.var.obj = value(coef.expr) - self.gurobi_model.ObjCon = value(self.constant.expr) - - gurobi_expr = None - for ndx, coef in enumerate(self.quadratic_coefs): - if value(coef.expr) != self.last_quadratic_coef_values[ndx]: - if gurobi_expr is None: - self.gurobi_model.update() - gurobi_expr = self.gurobi_model.getObjective() - current_coef_value = value(coef.expr) - incremental_coef_value = ( - current_coef_value - self.last_quadratic_coef_values[ndx] - ) - gurobi_expr += incremental_coef_value * coef.var1 * coef.var2 - self.last_quadratic_coef_values[ndx] = current_coef_value - return gurobi_expr - - -class _MutableQuadraticCoefficient: - def __init__(self): - self.expr = None - self.var1 = None - self.var2 = None - - -class GurobiPersistent( - GurobiSolverMixin, - PersistentSolverMixin, - PersistentSolverUtils, - PersistentSolverBase, -): - """ - Interface to Gurobi persistent - """ - - CONFIG = GurobiConfig() - _gurobipy_available = gurobipy_available - - def __init__(self, **kwds): - treat_fixed_vars_as_params = kwds.pop('treat_fixed_vars_as_params', True) - PersistentSolverBase.__init__(self, **kwds) - PersistentSolverUtils.__init__( - self, treat_fixed_vars_as_params=treat_fixed_vars_as_params - ) - self._register_env_client() - self._solver_model = None - self._symbol_map = SymbolMap() - self._labeler = None - self._pyomo_var_to_solver_var_map = {} - self._pyomo_con_to_solver_con_map = {} - self._solver_con_to_pyomo_con_map = {} - self._pyomo_sos_to_solver_sos_map = {} - self._range_constraints = OrderedSet() - self._mutable_helpers = {} - self._mutable_bounds = {} - self._mutable_quadratic_helpers = {} - self._mutable_objective = None - self._needs_updated = True - self._callback = None - self._callback_func = None - self._constraints_added_since_update = OrderedSet() - self._vars_added_since_update = ComponentSet() - self._last_results_object: Optional[Results] = None - - def release_license(self): - self._reinit() - self.__class__.release_license() - - def __del__(self): - if not python_is_shutting_down(): - self._release_env_client() - - @property - def symbol_map(self): - return self._symbol_map - - def _solve(self): - config = self._active_config - timer = config.timer - ostreams = [io.StringIO()] + config.tee - - with capture_output(TeeStream(*ostreams), capture_fd=False): - options = config.solver_options - - self._solver_model.setParam('LogToConsole', 1) - - if config.threads is not None: - self._solver_model.setParam('Threads', config.threads) - if config.time_limit is not None: - self._solver_model.setParam('TimeLimit', config.time_limit) - if config.rel_gap is not None: - self._solver_model.setParam('MIPGap', config.rel_gap) - if config.abs_gap is not None: - self._solver_model.setParam('MIPGapAbs', config.abs_gap) - - if config.use_mipstart: - for ( - pyomo_var_id, - gurobi_var, - ) in self._pyomo_var_to_solver_var_map.items(): - pyomo_var = self._vars[pyomo_var_id][0] - if pyomo_var.is_integer() and pyomo_var.value is not None: - self.set_var_attr(pyomo_var, 'Start', pyomo_var.value) - - for key, option in options.items(): - self._solver_model.setParam(key, option) - - timer.start('optimize') - self._solver_model.optimize(self._callback) - timer.stop('optimize') - - self._needs_updated = False - res = self._postsolve(timer) - res.solver_config = config - res.solver_name = 'Gurobi' - res.solver_version = self.version() - res.solver_log = ostreams[0].getvalue() - return res - - def _process_domain_and_bounds( - self, var, var_id, mutable_lbs, mutable_ubs, ndx, gurobipy_var - ): - _v, _lb, _ub, _fixed, _domain_interval, _value = self._vars[id(var)] - lb, ub, step = _domain_interval - if lb is None: - lb = -gurobipy.GRB.INFINITY - if ub is None: - ub = gurobipy.GRB.INFINITY - if step == 0: - vtype = gurobipy.GRB.CONTINUOUS - elif step == 1: - if lb == 0 and ub == 1: - vtype = gurobipy.GRB.BINARY - else: - vtype = gurobipy.GRB.INTEGER - else: - raise ValueError( - f'Unrecognized domain step: {step} (should be either 0 or 1)' - ) - if _fixed: - lb = _value - ub = _value - else: - if _lb is not None: - if not is_constant(_lb): - mutable_bound = _MutableLowerBound(NPV_MaxExpression((_lb, lb))) - if gurobipy_var is None: - mutable_lbs[ndx] = mutable_bound - else: - mutable_bound.var = gurobipy_var - self._mutable_bounds[var_id, 'lb'] = (var, mutable_bound) - lb = max(value(_lb), lb) - if _ub is not None: - if not is_constant(_ub): - mutable_bound = _MutableUpperBound(NPV_MinExpression((_ub, ub))) - if gurobipy_var is None: - mutable_ubs[ndx] = mutable_bound - else: - mutable_bound.var = gurobipy_var - self._mutable_bounds[var_id, 'ub'] = (var, mutable_bound) - ub = min(value(_ub), ub) - - return lb, ub, vtype - - def _add_variables(self, variables: List[VarData]): - var_names = [] - vtypes = [] - lbs = [] - ubs = [] - mutable_lbs = {} - mutable_ubs = {} - for ndx, var in enumerate(variables): - varname = self._symbol_map.getSymbol(var, self._labeler) - lb, ub, vtype = self._process_domain_and_bounds( - var, id(var), mutable_lbs, mutable_ubs, ndx, None - ) - var_names.append(varname) - vtypes.append(vtype) - lbs.append(lb) - ubs.append(ub) - - gurobi_vars = self._solver_model.addVars( - len(variables), lb=lbs, ub=ubs, vtype=vtypes, name=var_names - ) - - for ndx, pyomo_var in enumerate(variables): - gurobi_var = gurobi_vars[ndx] - self._pyomo_var_to_solver_var_map[id(pyomo_var)] = gurobi_var - for ndx, mutable_bound in mutable_lbs.items(): - mutable_bound.var = gurobi_vars[ndx] - for ndx, mutable_bound in mutable_ubs.items(): - mutable_bound.var = gurobi_vars[ndx] - self._vars_added_since_update.update(variables) - self._needs_updated = True - - def _add_parameters(self, params: List[ParamData]): - pass - - def _reinit(self): - saved_config = self.config - saved_tmp_config = self._active_config - self.__init__(treat_fixed_vars_as_params=self._treat_fixed_vars_as_params) - # Note that __init__ registers a new env client, so we need to - # release it here: - self._release_env_client() - self.config = saved_config - self._active_config = saved_tmp_config - - def set_instance(self, model): - if self._last_results_object is not None: - self._last_results_object.solution_loader.invalidate() - if not self.available(): - c = self.__class__ - raise ApplicationError( - f'Solver {c.__module__}.{c.__qualname__} is not available ' - f'({self.available()}).' - ) - self._reinit() - self._model = model - - if self.config.symbolic_solver_labels: - self._labeler = TextLabeler() - else: - self._labeler = NumericLabeler('x') - - self._solver_model = gurobipy.Model(name=model.name or '', env=self.env()) - - self.add_block(model) - if self._objective is None: - self.set_objective(None) - - def _get_expr_from_pyomo_expr(self, expr): - mutable_linear_coefficients = [] - mutable_quadratic_coefficients = [] - repn = generate_standard_repn(expr, quadratic=True, compute_values=False) - - degree = repn.polynomial_degree() - if (degree is None) or (degree > 2): - raise IncompatibleModelError( - f'GurobiAuto does not support expressions of degree {degree}.' - ) - - if len(repn.linear_vars) > 0: - linear_coef_vals = [] - for ndx, coef in enumerate(repn.linear_coefs): - if not is_constant(coef): - mutable_linear_coefficient = _MutableLinearCoefficient() - mutable_linear_coefficient.expr = coef - mutable_linear_coefficient.var = self._pyomo_var_to_solver_var_map[ - id(repn.linear_vars[ndx]) - ] - mutable_linear_coefficients.append(mutable_linear_coefficient) - linear_coef_vals.append(value(coef)) - new_expr = gurobipy.LinExpr( - linear_coef_vals, - [self._pyomo_var_to_solver_var_map[id(i)] for i in repn.linear_vars], - ) - else: - new_expr = 0.0 - - for ndx, v in enumerate(repn.quadratic_vars): - x, y = v - gurobi_x = self._pyomo_var_to_solver_var_map[id(x)] - gurobi_y = self._pyomo_var_to_solver_var_map[id(y)] - coef = repn.quadratic_coefs[ndx] - if not is_constant(coef): - mutable_quadratic_coefficient = _MutableQuadraticCoefficient() - mutable_quadratic_coefficient.expr = coef - mutable_quadratic_coefficient.var1 = gurobi_x - mutable_quadratic_coefficient.var2 = gurobi_y - mutable_quadratic_coefficients.append(mutable_quadratic_coefficient) - coef_val = value(coef) - new_expr += coef_val * gurobi_x * gurobi_y - - return ( - new_expr, - repn.constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) - - def _add_constraints(self, cons: List[ConstraintData]): - for con in cons: - conname = self._symbol_map.getSymbol(con, self._labeler) - ( - gurobi_expr, - repn_constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) = self._get_expr_from_pyomo_expr(con.body) - - if ( - gurobi_expr.__class__ in {gurobipy.LinExpr, gurobipy.Var} - or gurobi_expr.__class__ in native_numeric_types - ): - if con.equality: - rhs_expr = con.lower - repn_constant - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addLConstr( - gurobi_expr, gurobipy.GRB.EQUAL, rhs_val, name=conname - ) - if not is_constant(rhs_expr): - mutable_constant = _MutableConstant() - mutable_constant.expr = rhs_expr - mutable_constant.con = gurobipy_con - self._mutable_helpers[con] = [mutable_constant] - elif con.has_lb() and con.has_ub(): - lhs_expr = con.lower - repn_constant - rhs_expr = con.upper - repn_constant - lhs_val = value(lhs_expr) - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addRange( - gurobi_expr, lhs_val, rhs_val, name=conname - ) - self._range_constraints.add(con) - if not is_constant(lhs_expr) or not is_constant(rhs_expr): - mutable_range_constant = _MutableRangeConstant() - mutable_range_constant.lhs_expr = lhs_expr - mutable_range_constant.rhs_expr = rhs_expr - mutable_range_constant.con = gurobipy_con - mutable_range_constant.slack_name = 'Rg' + conname - mutable_range_constant.gurobi_model = self._solver_model - self._mutable_helpers[con] = [mutable_range_constant] - elif con.has_lb(): - rhs_expr = con.lower - repn_constant - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addLConstr( - gurobi_expr, gurobipy.GRB.GREATER_EQUAL, rhs_val, name=conname - ) - if not is_constant(rhs_expr): - mutable_constant = _MutableConstant() - mutable_constant.expr = rhs_expr - mutable_constant.con = gurobipy_con - self._mutable_helpers[con] = [mutable_constant] - elif con.has_ub(): - rhs_expr = con.upper - repn_constant - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addLConstr( - gurobi_expr, gurobipy.GRB.LESS_EQUAL, rhs_val, name=conname - ) - if not is_constant(rhs_expr): - mutable_constant = _MutableConstant() - mutable_constant.expr = rhs_expr - mutable_constant.con = gurobipy_con - self._mutable_helpers[con] = [mutable_constant] - else: - raise ValueError( - "Constraint does not have a lower " - f"or an upper bound: {con} \n" - ) - for tmp in mutable_linear_coefficients: - tmp.con = gurobipy_con - tmp.gurobi_model = self._solver_model - if len(mutable_linear_coefficients) > 0: - if con not in self._mutable_helpers: - self._mutable_helpers[con] = mutable_linear_coefficients - else: - self._mutable_helpers[con].extend(mutable_linear_coefficients) - elif gurobi_expr.__class__ is gurobipy.QuadExpr: - if con.equality: - rhs_expr = con.lower - repn_constant - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addQConstr( - gurobi_expr, gurobipy.GRB.EQUAL, rhs_val, name=conname - ) - elif con.has_lb() and con.has_ub(): - raise NotImplementedError( - 'Quadratic range constraints are not supported' - ) - elif con.has_lb(): - rhs_expr = con.lower - repn_constant - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addQConstr( - gurobi_expr, gurobipy.GRB.GREATER_EQUAL, rhs_val, name=conname - ) - elif con.has_ub(): - rhs_expr = con.upper - repn_constant - rhs_val = value(rhs_expr) - gurobipy_con = self._solver_model.addQConstr( - gurobi_expr, gurobipy.GRB.LESS_EQUAL, rhs_val, name=conname - ) - else: - raise ValueError( - "Constraint does not have a lower " - f"or an upper bound: {con} \n" - ) - if ( - len(mutable_linear_coefficients) > 0 - or len(mutable_quadratic_coefficients) > 0 - or not is_constant(repn_constant) - ): - mutable_constant = _MutableConstant() - mutable_constant.expr = rhs_expr - mutable_quadratic_constraint = _MutableQuadraticConstraint( - self._solver_model, - gurobipy_con, - mutable_constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) - self._mutable_quadratic_helpers[con] = mutable_quadratic_constraint - else: - raise ValueError( - f'Unrecognized Gurobi expression type: {str(gurobi_expr.__class__)}' - ) - - self._pyomo_con_to_solver_con_map[con] = gurobipy_con - self._solver_con_to_pyomo_con_map[id(gurobipy_con)] = con - self._constraints_added_since_update.update(cons) - self._needs_updated = True - - def _add_sos_constraints(self, cons: List[SOSConstraintData]): - for con in cons: - conname = self._symbol_map.getSymbol(con, self._labeler) - level = con.level - if level == 1: - sos_type = gurobipy.GRB.SOS_TYPE1 - elif level == 2: - sos_type = gurobipy.GRB.SOS_TYPE2 - else: - raise ValueError( - f"Solver does not support SOS level {level} constraints" - ) - - gurobi_vars = [] - weights = [] - - for v, w in con.get_items(): - v_id = id(v) - gurobi_vars.append(self._pyomo_var_to_solver_var_map[v_id]) - weights.append(w) - - gurobipy_con = self._solver_model.addSOS(sos_type, gurobi_vars, weights) - self._pyomo_sos_to_solver_sos_map[con] = gurobipy_con - self._constraints_added_since_update.update(cons) - self._needs_updated = True - - def _remove_constraints(self, cons: List[ConstraintData]): - for con in cons: - if con in self._constraints_added_since_update: - self._update_gurobi_model() - solver_con = self._pyomo_con_to_solver_con_map[con] - self._solver_model.remove(solver_con) - self._symbol_map.removeSymbol(con) - del self._pyomo_con_to_solver_con_map[con] - del self._solver_con_to_pyomo_con_map[id(solver_con)] - self._range_constraints.discard(con) - self._mutable_helpers.pop(con, None) - self._mutable_quadratic_helpers.pop(con, None) - self._needs_updated = True - - def _remove_sos_constraints(self, cons: List[SOSConstraintData]): - for con in cons: - if con in self._constraints_added_since_update: - self._update_gurobi_model() - solver_sos_con = self._pyomo_sos_to_solver_sos_map[con] - self._solver_model.remove(solver_sos_con) - self._symbol_map.removeSymbol(con) - del self._pyomo_sos_to_solver_sos_map[con] - self._needs_updated = True - - def _remove_variables(self, variables: List[VarData]): - for var in variables: - v_id = id(var) - if var in self._vars_added_since_update: - self._update_gurobi_model() - solver_var = self._pyomo_var_to_solver_var_map[v_id] - self._solver_model.remove(solver_var) - self._symbol_map.removeSymbol(var) - del self._pyomo_var_to_solver_var_map[v_id] - self._mutable_bounds.pop(v_id, None) - self._needs_updated = True - - def _remove_parameters(self, params: List[ParamData]): - pass - - def _update_variables(self, variables: List[VarData]): - for var in variables: - var_id = id(var) - if var_id not in self._pyomo_var_to_solver_var_map: - raise ValueError( - f'The Var provided to update_var needs to be added first: {var}' - ) - self._mutable_bounds.pop((var_id, 'lb'), None) - self._mutable_bounds.pop((var_id, 'ub'), None) - gurobipy_var = self._pyomo_var_to_solver_var_map[var_id] - lb, ub, vtype = self._process_domain_and_bounds( - var, var_id, None, None, None, gurobipy_var - ) - gurobipy_var.setAttr('lb', lb) - gurobipy_var.setAttr('ub', ub) - gurobipy_var.setAttr('vtype', vtype) - self._needs_updated = True - - def update_parameters(self): - for con, helpers in self._mutable_helpers.items(): - for helper in helpers: - helper.update() - for k, (v, helper) in self._mutable_bounds.items(): - helper.update() - - for con, helper in self._mutable_quadratic_helpers.items(): - if con in self._constraints_added_since_update: - self._update_gurobi_model() - gurobi_con = helper.con - new_gurobi_expr = helper.get_updated_expression() - new_rhs = helper.get_updated_rhs() - new_sense = gurobi_con.qcsense - pyomo_con = self._solver_con_to_pyomo_con_map[id(gurobi_con)] - name = self._symbol_map.getSymbol(pyomo_con, self._labeler) - self._solver_model.remove(gurobi_con) - new_con = self._solver_model.addQConstr( - new_gurobi_expr, new_sense, new_rhs, name=name - ) - self._pyomo_con_to_solver_con_map[id(pyomo_con)] = new_con - del self._solver_con_to_pyomo_con_map[id(gurobi_con)] - self._solver_con_to_pyomo_con_map[id(new_con)] = pyomo_con - helper.con = new_con - self._constraints_added_since_update.add(con) - - helper = self._mutable_objective - pyomo_obj = self._objective - new_gurobi_expr = helper.get_updated_expression() - if new_gurobi_expr is not None: - if pyomo_obj.sense == minimize: - sense = gurobipy.GRB.MINIMIZE - else: - sense = gurobipy.GRB.MAXIMIZE - self._solver_model.setObjective(new_gurobi_expr, sense=sense) - - def _set_objective(self, obj): - if obj is None: - sense = gurobipy.GRB.MINIMIZE - gurobi_expr = 0 - repn_constant = 0 - mutable_linear_coefficients = [] - mutable_quadratic_coefficients = [] - else: - if obj.sense == minimize: - sense = gurobipy.GRB.MINIMIZE - elif obj.sense == maximize: - sense = gurobipy.GRB.MAXIMIZE - else: - raise ValueError(f'Objective sense is not recognized: {obj.sense}') - - ( - gurobi_expr, - repn_constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) = self._get_expr_from_pyomo_expr(obj.expr) - - mutable_constant = _MutableConstant() - mutable_constant.expr = repn_constant - mutable_objective = _MutableObjective( - self._solver_model, - mutable_constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) - self._mutable_objective = mutable_objective - - # These two lines are needed as a workaround - # see PR #2454 - self._solver_model.setObjective(0) - self._solver_model.update() - - self._solver_model.setObjective(gurobi_expr + value(repn_constant), sense=sense) - self._needs_updated = True - - def _postsolve(self, timer: HierarchicalTimer): - config = self._active_config - - gprob = self._solver_model - grb = gurobipy.GRB - status = gprob.Status - - results = Results() - results.solution_loader = GurobiSolutionLoader(self) - results.timing_info.gurobi_time = gprob.Runtime - - if gprob.SolCount > 0: - if status == grb.OPTIMAL: - results.solution_status = SolutionStatus.optimal - else: - results.solution_status = SolutionStatus.feasible - else: - results.solution_status = SolutionStatus.noSolution - - if status == grb.LOADED: # problem is loaded, but no solution - results.termination_condition = TerminationCondition.unknown - elif status == grb.OPTIMAL: # optimal - results.termination_condition = ( - TerminationCondition.convergenceCriteriaSatisfied - ) - elif status == grb.INFEASIBLE: - results.termination_condition = TerminationCondition.provenInfeasible - elif status == grb.INF_OR_UNBD: - results.termination_condition = TerminationCondition.infeasibleOrUnbounded - elif status == grb.UNBOUNDED: - results.termination_condition = TerminationCondition.unbounded - elif status == grb.CUTOFF: - results.termination_condition = TerminationCondition.objectiveLimit - elif status == grb.ITERATION_LIMIT: - results.termination_condition = TerminationCondition.iterationLimit - elif status == grb.NODE_LIMIT: - results.termination_condition = TerminationCondition.iterationLimit - elif status == grb.TIME_LIMIT: - results.termination_condition = TerminationCondition.maxTimeLimit - elif status == grb.SOLUTION_LIMIT: - results.termination_condition = TerminationCondition.unknown - elif status == grb.INTERRUPTED: - results.termination_condition = TerminationCondition.interrupted - elif status == grb.NUMERIC: - results.termination_condition = TerminationCondition.unknown - elif status == grb.SUBOPTIMAL: - results.termination_condition = TerminationCondition.unknown - elif status == grb.USER_OBJ_LIMIT: - results.termination_condition = TerminationCondition.objectiveLimit - else: - results.termination_condition = TerminationCondition.unknown - - if ( - results.termination_condition - != TerminationCondition.convergenceCriteriaSatisfied - and config.raise_exception_on_nonoptimal_result - ): - raise NoOptimalSolutionError() - - results.incumbent_objective = None - results.objective_bound = None - if self._objective is not None: - try: - results.incumbent_objective = gprob.ObjVal - except (gurobipy.GurobiError, AttributeError): - results.incumbent_objective = None - try: - results.objective_bound = gprob.ObjBound - except (gurobipy.GurobiError, AttributeError): - if self._objective.sense == minimize: - results.objective_bound = -math.inf - else: - results.objective_bound = math.inf - - if results.incumbent_objective is not None and not math.isfinite( - results.incumbent_objective - ): - results.incumbent_objective = None - - results.iteration_count = gprob.getAttr('IterCount') - - timer.start('load solution') - if config.load_solutions: - if gprob.SolCount > 0: - self._load_vars() - else: - raise NoFeasibleSolutionError() - timer.stop('load solution') - - return results - - def _load_suboptimal_mip_solution(self, vars_to_load, solution_number): - if ( - self.get_model_attr('NumIntVars') == 0 - and self.get_model_attr('NumBinVars') == 0 - ): - raise ValueError( - 'Cannot obtain suboptimal solutions for a continuous model' - ) - var_map = self._pyomo_var_to_solver_var_map - ref_vars = self._referenced_variables - original_solution_number = self.get_gurobi_param_info('SolutionNumber')[2] - self.set_gurobi_param('SolutionNumber', solution_number) - gurobi_vars_to_load = [var_map[pyomo_var] for pyomo_var in vars_to_load] - vals = self._solver_model.getAttr("Xn", gurobi_vars_to_load) - res = ComponentMap() - for var_id, val in zip(vars_to_load, vals): - using_cons, using_sos, using_obj = ref_vars[var_id] - if using_cons or using_sos or (using_obj is not None): - res[self._vars[var_id][0]] = val - self.set_gurobi_param('SolutionNumber', original_solution_number) - return res - - def _load_vars(self, vars_to_load=None, solution_number=0): - for v, val in self._get_primals( - vars_to_load=vars_to_load, solution_number=solution_number - ).items(): - v.set_value(val, skip_validation=True) - StaleFlagManager.mark_all_as_stale(delayed=True) - - def _get_primals(self, vars_to_load=None, solution_number=0): - if self._needs_updated: - self._update_gurobi_model() # this is needed to ensure that solutions cannot be loaded after the model has been changed - - if self._solver_model.SolCount == 0: - raise NoSolutionError() - - var_map = self._pyomo_var_to_solver_var_map - ref_vars = self._referenced_variables - if vars_to_load is None: - vars_to_load = self._pyomo_var_to_solver_var_map.keys() - else: - vars_to_load = [id(v) for v in vars_to_load] - - if solution_number != 0: - return self._load_suboptimal_mip_solution( - vars_to_load=vars_to_load, solution_number=solution_number - ) - - gurobi_vars_to_load = [var_map[pyomo_var_id] for pyomo_var_id in vars_to_load] - vals = self._solver_model.getAttr("X", gurobi_vars_to_load) - - res = ComponentMap() - for var_id, val in zip(vars_to_load, vals): - using_cons, using_sos, using_obj = ref_vars[var_id] - if using_cons or using_sos or (using_obj is not None): - res[self._vars[var_id][0]] = val - return res - - def _get_reduced_costs(self, vars_to_load=None): - if self._needs_updated: - self._update_gurobi_model() - - if self._solver_model.Status != gurobipy.GRB.OPTIMAL: - raise NoReducedCostsError() - - var_map = self._pyomo_var_to_solver_var_map - ref_vars = self._referenced_variables - res = ComponentMap() - if vars_to_load is None: - vars_to_load = self._pyomo_var_to_solver_var_map.keys() - else: - vars_to_load = [id(v) for v in vars_to_load] - - gurobi_vars_to_load = [var_map[pyomo_var_id] for pyomo_var_id in vars_to_load] - vals = self._solver_model.getAttr("Rc", gurobi_vars_to_load) - - for var_id, val in zip(vars_to_load, vals): - using_cons, using_sos, using_obj = ref_vars[var_id] - if using_cons or using_sos or (using_obj is not None): - res[self._vars[var_id][0]] = val - - return res - - def _get_duals(self, cons_to_load=None): - if self._needs_updated: - self._update_gurobi_model() - - if self._solver_model.Status != gurobipy.GRB.OPTIMAL: - raise NoDualsError() - - con_map = self._pyomo_con_to_solver_con_map - reverse_con_map = self._solver_con_to_pyomo_con_map - dual = {} - - if cons_to_load is None: - linear_cons_to_load = self._solver_model.getConstrs() - quadratic_cons_to_load = self._solver_model.getQConstrs() - else: - gurobi_cons_to_load = OrderedSet( - [con_map[pyomo_con] for pyomo_con in cons_to_load] - ) - linear_cons_to_load = list( - gurobi_cons_to_load.intersection( - OrderedSet(self._solver_model.getConstrs()) - ) - ) - quadratic_cons_to_load = list( - gurobi_cons_to_load.intersection( - OrderedSet(self._solver_model.getQConstrs()) - ) - ) - linear_vals = self._solver_model.getAttr("Pi", linear_cons_to_load) - quadratic_vals = self._solver_model.getAttr("QCPi", quadratic_cons_to_load) - - for gurobi_con, val in zip(linear_cons_to_load, linear_vals): - pyomo_con = reverse_con_map[id(gurobi_con)] - dual[pyomo_con] = val - for gurobi_con, val in zip(quadratic_cons_to_load, quadratic_vals): - pyomo_con = reverse_con_map[id(gurobi_con)] - dual[pyomo_con] = val - - return dual - - def update(self, timer: HierarchicalTimer = None): - if self._needs_updated: - self._update_gurobi_model() - super().update(timer=timer) - self._update_gurobi_model() - - def _update_gurobi_model(self): - self._solver_model.update() - self._constraints_added_since_update = OrderedSet() - self._vars_added_since_update = ComponentSet() - self._needs_updated = False - - def get_model_attr(self, attr): - """ - Get the value of an attribute on the Gurobi model. - - Parameters - ---------- - attr: str - The attribute to get. See Gurobi documentation for descriptions of the attributes. - """ - if self._needs_updated: - self._update_gurobi_model() - return self._solver_model.getAttr(attr) - - def write(self, filename): - """ - Write the model to a file (e.g., and lp file). - - Parameters - ---------- - filename: str - Name of the file to which the model should be written. - """ - self._solver_model.write(filename) - self._constraints_added_since_update = OrderedSet() - self._vars_added_since_update = ComponentSet() - self._needs_updated = False - - def set_linear_constraint_attr(self, con, attr, val): - """ - Set the value of an attribute on a gurobi linear constraint. - - Parameters - ---------- - con: pyomo.core.base.constraint.ConstraintData - The pyomo constraint for which the corresponding gurobi constraint attribute - should be modified. - attr: str - The attribute to be modified. Options are: - CBasis - DStart - Lazy - val: any - See gurobi documentation for acceptable values. - """ - if attr in {'Sense', 'RHS', 'ConstrName'}: - raise ValueError( - f'Linear constraint attr {attr} cannot be set with' - ' the set_linear_constraint_attr method. Please use' - ' the remove_constraint and add_constraint methods.' - ) - self._pyomo_con_to_solver_con_map[con].setAttr(attr, val) - self._needs_updated = True - - def set_var_attr(self, var, attr, val): - """ - Set the value of an attribute on a gurobi variable. - - Parameters - ---------- - var: pyomo.core.base.var.VarData - The pyomo var for which the corresponding gurobi var attribute - should be modified. - attr: str - The attribute to be modified. Options are: - Start - VarHintVal - VarHintPri - BranchPriority - VBasis - PStart - val: any - See gurobi documentation for acceptable values. - """ - if attr in {'LB', 'UB', 'VType', 'VarName'}: - raise ValueError( - f'Var attr {attr} cannot be set with' - ' the set_var_attr method. Please use' - ' the update_var method.' - ) - if attr == 'Obj': - raise ValueError( - 'Var attr Obj cannot be set with' - ' the set_var_attr method. Please use' - ' the set_objective method.' - ) - self._pyomo_var_to_solver_var_map[id(var)].setAttr(attr, val) - self._needs_updated = True - - def get_var_attr(self, var, attr): - """ - Get the value of an attribute on a gurobi var. - - Parameters - ---------- - var: pyomo.core.base.var.VarData - The pyomo var for which the corresponding gurobi var attribute - should be retrieved. - attr: str - The attribute to get. See gurobi documentation - """ - if self._needs_updated: - self._update_gurobi_model() - return self._pyomo_var_to_solver_var_map[id(var)].getAttr(attr) - - def get_linear_constraint_attr(self, con, attr): - """ - Get the value of an attribute on a gurobi linear constraint. - - Parameters - ---------- - con: pyomo.core.base.constraint.ConstraintData - The pyomo constraint for which the corresponding gurobi constraint attribute - should be retrieved. - attr: str - The attribute to get. See the Gurobi documentation - """ - if self._needs_updated: - self._update_gurobi_model() - return self._pyomo_con_to_solver_con_map[con].getAttr(attr) - - def get_sos_attr(self, con, attr): - """ - Get the value of an attribute on a gurobi sos constraint. - - Parameters - ---------- - con: pyomo.core.base.sos.SOSConstraintData - The pyomo SOS constraint for which the corresponding gurobi SOS constraint attribute - should be retrieved. - attr: str - The attribute to get. See the Gurobi documentation - """ - if self._needs_updated: - self._update_gurobi_model() - return self._pyomo_sos_to_solver_sos_map[con].getAttr(attr) - - def get_quadratic_constraint_attr(self, con, attr): - """ - Get the value of an attribute on a gurobi quadratic constraint. - - Parameters - ---------- - con: pyomo.core.base.constraint.ConstraintData - The pyomo constraint for which the corresponding gurobi constraint attribute - should be retrieved. - attr: str - The attribute to get. See the Gurobi documentation - """ - if self._needs_updated: - self._update_gurobi_model() - return self._pyomo_con_to_solver_con_map[con].getAttr(attr) - - def set_gurobi_param(self, param, val): - """ - Set a gurobi parameter. - - Parameters - ---------- - param: str - The gurobi parameter to set. Options include any gurobi parameter. - Please see the Gurobi documentation for options. - val: any - The value to set the parameter to. See Gurobi documentation for possible values. - """ - self._solver_model.setParam(param, val) - - def get_gurobi_param_info(self, param): - """ - Get information about a gurobi parameter. - - Parameters - ---------- - param: str - The gurobi parameter to get info for. See Gurobi documentation for possible options. - - Returns - ------- - six-tuple containing the parameter name, type, value, minimum value, maximum value, and default value. - """ - return self._solver_model.getParamInfo(param) - - def _intermediate_callback(self): - def f(gurobi_model, where): - self._callback_func(self._model, self, where) - - return f - - def set_callback(self, func=None): - """ - Specify a callback for gurobi to use. - - Parameters - ---------- - func: function - The function to call. The function should have three arguments. The first will be the pyomo model being - solved. The second will be the GurobiPersistent instance. The third will be an enum member of - gurobipy.GRB.Callback. This will indicate where in the branch and bound algorithm gurobi is at. For - example, suppose we want to solve - - .. math:: - - min 2*x + y - - s.t. - - y >= (x-2)**2 - - 0 <= x <= 4 - - y >= 0 - - y integer - - as an MILP using extended cutting planes in callbacks. - - >>> from gurobipy import GRB # doctest:+SKIP - >>> import pyomo.environ as pyo - >>> from pyomo.core.expr.taylor_series import taylor_series_expansion - >>> from pyomo.contrib import appsi - >>> - >>> m = pyo.ConcreteModel() - >>> m.x = pyo.Var(bounds=(0, 4)) - >>> m.y = pyo.Var(within=pyo.Integers, bounds=(0, None)) - >>> m.obj = pyo.Objective(expr=2*m.x + m.y) - >>> m.cons = pyo.ConstraintList() # for the cutting planes - >>> - >>> def _add_cut(xval): - ... # a function to generate the cut - ... m.x.value = xval - ... return m.cons.add(m.y >= taylor_series_expansion((m.x - 2)**2)) - ... - >>> _c = _add_cut(0) # start with 2 cuts at the bounds of x - >>> _c = _add_cut(4) # this is an arbitrary choice - >>> - >>> opt = appsi.solvers.Gurobi() - >>> opt.config.stream_solver = True - >>> opt.set_instance(m) # doctest:+SKIP - >>> opt.gurobi_options['PreCrush'] = 1 - >>> opt.gurobi_options['LazyConstraints'] = 1 - >>> - >>> def my_callback(cb_m, cb_opt, cb_where): - ... if cb_where == GRB.Callback.MIPSOL: - ... cb_opt.cbGetSolution(variables=[m.x, m.y]) - ... if m.y.value < (m.x.value - 2)**2 - 1e-6: - ... cb_opt.cbLazy(_add_cut(m.x.value)) - ... - >>> opt.set_callback(my_callback) - >>> res = opt.solve(m) # doctest:+SKIP - - """ - if func is not None: - self._callback_func = func - self._callback = self._intermediate_callback() - else: - self._callback = None - self._callback_func = None - - def cbCut(self, con): - """ - Add a cut within a callback. - - Parameters - ---------- - con: pyomo.core.base.constraint.ConstraintData - The cut to add - """ - if not con.active: - raise ValueError('cbCut expected an active constraint.') - - if is_fixed(con.body): - raise ValueError('cbCut expected a non-trivial constraint') - - ( - gurobi_expr, - repn_constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) = self._get_expr_from_pyomo_expr(con.body) - - if con.has_lb(): - if con.has_ub(): - raise ValueError('Range constraints are not supported in cbCut.') - if not is_fixed(con.lower): - raise ValueError(f'Lower bound of constraint {con} is not constant.') - if con.has_ub(): - if not is_fixed(con.upper): - raise ValueError(f'Upper bound of constraint {con} is not constant.') - - if con.equality: - self._solver_model.cbCut( - lhs=gurobi_expr, - sense=gurobipy.GRB.EQUAL, - rhs=value(con.lower - repn_constant), - ) - elif con.has_lb() and (value(con.lower) > -float('inf')): - self._solver_model.cbCut( - lhs=gurobi_expr, - sense=gurobipy.GRB.GREATER_EQUAL, - rhs=value(con.lower - repn_constant), - ) - elif con.has_ub() and (value(con.upper) < float('inf')): - self._solver_model.cbCut( - lhs=gurobi_expr, - sense=gurobipy.GRB.LESS_EQUAL, - rhs=value(con.upper - repn_constant), - ) - else: - raise ValueError( - f'Constraint does not have a lower or an upper bound {con} \n' - ) - - def cbGet(self, what): - return self._solver_model.cbGet(what) - - def cbGetNodeRel(self, variables): - """ - Parameters - ---------- - variables: Var or iterable of Var - """ - if not isinstance(variables, Iterable): - variables = [variables] - gurobi_vars = [self._pyomo_var_to_solver_var_map[id(i)] for i in variables] - var_values = self._solver_model.cbGetNodeRel(gurobi_vars) - for i, v in enumerate(variables): - v.set_value(var_values[i], skip_validation=True) - - def cbGetSolution(self, variables): - """ - Parameters - ---------- - variables: iterable of vars - """ - if not isinstance(variables, Iterable): - variables = [variables] - gurobi_vars = [self._pyomo_var_to_solver_var_map[id(i)] for i in variables] - var_values = self._solver_model.cbGetSolution(gurobi_vars) - for i, v in enumerate(variables): - v.set_value(var_values[i], skip_validation=True) - - def cbLazy(self, con): - """ - Parameters - ---------- - con: pyomo.core.base.constraint.ConstraintData - The lazy constraint to add - """ - if not con.active: - raise ValueError('cbLazy expected an active constraint.') - - if is_fixed(con.body): - raise ValueError('cbLazy expected a non-trivial constraint') - - ( - gurobi_expr, - repn_constant, - mutable_linear_coefficients, - mutable_quadratic_coefficients, - ) = self._get_expr_from_pyomo_expr(con.body) - - if con.has_lb(): - if con.has_ub(): - raise ValueError('Range constraints are not supported in cbLazy.') - if not is_fixed(con.lower): - raise ValueError(f'Lower bound of constraint {con} is not constant.') - if con.has_ub(): - if not is_fixed(con.upper): - raise ValueError(f'Upper bound of constraint {con} is not constant.') - - if con.equality: - self._solver_model.cbLazy( - lhs=gurobi_expr, - sense=gurobipy.GRB.EQUAL, - rhs=value(con.lower - repn_constant), - ) - elif con.has_lb() and (value(con.lower) > -float('inf')): - self._solver_model.cbLazy( - lhs=gurobi_expr, - sense=gurobipy.GRB.GREATER_EQUAL, - rhs=value(con.lower - repn_constant), - ) - elif con.has_ub() and (value(con.upper) < float('inf')): - self._solver_model.cbLazy( - lhs=gurobi_expr, - sense=gurobipy.GRB.LESS_EQUAL, - rhs=value(con.upper - repn_constant), - ) - else: - raise ValueError( - f'Constraint does not have a lower or an upper bound {con} \n' - ) - - def cbSetSolution(self, variables, solution): - if not isinstance(variables, Iterable): - variables = [variables] - gurobi_vars = [self._pyomo_var_to_solver_var_map[id(i)] for i in variables] - self._solver_model.cbSetSolution(gurobi_vars, solution) - - def cbUseSolution(self): - return self._solver_model.cbUseSolution() - - def reset(self): - self._solver_model.reset() diff --git a/pyomo/contrib/solver/tests/solvers/test_gurobi_persistent.py b/pyomo/contrib/solver/tests/solvers/test_gurobi_persistent.py index 8703ae9edff..96cd1498956 100644 --- a/pyomo/contrib/solver/tests/solvers/test_gurobi_persistent.py +++ b/pyomo/contrib/solver/tests/solvers/test_gurobi_persistent.py @@ -11,7 +11,7 @@ import pyomo.common.unittest as unittest import pyomo.environ as pyo -from pyomo.contrib.solver.solvers.gurobi_persistent import GurobiPersistent +from pyomo.contrib.solver.solvers.gurobi.gurobi_persistent import GurobiPersistent from pyomo.contrib.solver.common.results import SolutionStatus from pyomo.core.expr.taylor_series import taylor_series_expansion @@ -471,11 +471,11 @@ def test_solution_number(self): res = opt.solve(m) num_solutions = opt.get_model_attr('SolCount') self.assertEqual(num_solutions, 3) - res.solution_loader.load_vars(solution_number=0) + res.solution_loader.load_vars(solution_id=0) self.assertAlmostEqual(pyo.value(m.obj.expr), 6.431184939357673) - res.solution_loader.load_vars(solution_number=1) + res.solution_loader.load_vars(solution_id=1) self.assertAlmostEqual(pyo.value(m.obj.expr), 6.584793218502477) - res.solution_loader.load_vars(solution_number=2) + res.solution_loader.load_vars(solution_id=2) self.assertAlmostEqual(pyo.value(m.obj.expr), 6.592304628123309) def test_zero_time_limit(self): @@ -496,16 +496,14 @@ def test_zero_time_limit(self): self.assertIsNone(res.incumbent_objective) -class TestManualModel(unittest.TestCase): +class TestManualMode(unittest.TestCase): def setUp(self): opt = GurobiPersistent() - opt.config.auto_updates.check_for_new_or_removed_params = False - opt.config.auto_updates.check_for_new_or_removed_vars = False - opt.config.auto_updates.check_for_new_or_removed_constraints = False - opt.config.auto_updates.update_parameters = False - opt.config.auto_updates.update_vars = False - opt.config.auto_updates.update_constraints = False - opt.config.auto_updates.update_named_expressions = False + opt.auto_updates.check_for_new_or_removed_constraints = False + opt.auto_updates.update_parameters = False + opt.auto_updates.update_vars = False + opt.auto_updates.update_constraints = False + opt.auto_updates.update_named_expressions = False self.opt = opt def test_basics(self): @@ -603,16 +601,13 @@ def test_update1(self): opt = self.opt opt.set_instance(m) - self.assertEqual(opt._solver_model.getAttr('NumQConstrs'), 1) + self.assertEqual(opt.get_model_attr('NumQConstrs'), 1) opt.remove_constraints([m.c1]) - opt.update() - self.assertEqual(opt._solver_model.getAttr('NumQConstrs'), 0) + self.assertEqual(opt.get_model_attr('NumQConstrs'), 0) opt.add_constraints([m.c1]) - self.assertEqual(opt._solver_model.getAttr('NumQConstrs'), 0) - opt.update() - self.assertEqual(opt._solver_model.getAttr('NumQConstrs'), 1) + self.assertEqual(opt.get_model_attr('NumQConstrs'), 1) def test_update2(self): m = pyo.ConcreteModel() @@ -625,16 +620,13 @@ def test_update2(self): opt = self.opt opt.config.symbolic_solver_labels = True opt.set_instance(m) - self.assertEqual(opt._solver_model.getAttr('NumConstrs'), 1) + self.assertEqual(opt.get_model_attr('NumConstrs'), 1) opt.remove_constraints([m.c2]) - opt.update() - self.assertEqual(opt._solver_model.getAttr('NumConstrs'), 0) + self.assertEqual(opt.get_model_attr('NumConstrs'), 0) opt.add_constraints([m.c2]) - self.assertEqual(opt._solver_model.getAttr('NumConstrs'), 0) - opt.update() - self.assertEqual(opt._solver_model.getAttr('NumConstrs'), 1) + self.assertEqual(opt.get_model_attr('NumConstrs'), 1) def test_update3(self): m = pyo.ConcreteModel() @@ -684,16 +676,13 @@ def test_update5(self): opt = self.opt opt.set_instance(m) - self.assertEqual(opt._solver_model.getAttr('NumSOS'), 1) + self.assertEqual(opt.get_model_attr('NumSOS'), 1) opt.remove_sos_constraints([m.c1]) - opt.update() - self.assertEqual(opt._solver_model.getAttr('NumSOS'), 0) + self.assertEqual(opt.get_model_attr('NumSOS'), 0) opt.add_sos_constraints([m.c1]) - self.assertEqual(opt._solver_model.getAttr('NumSOS'), 0) - opt.update() - self.assertEqual(opt._solver_model.getAttr('NumSOS'), 1) + self.assertEqual(opt.get_model_attr('NumSOS'), 1) def test_update6(self): m = pyo.ConcreteModel() diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 5ab36554061..6965147d167 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -30,8 +30,11 @@ from pyomo.contrib.solver.common.base import SolverBase from pyomo.contrib.solver.common.factory import SolverFactory from pyomo.contrib.solver.solvers.ipopt import Ipopt -from pyomo.contrib.solver.solvers.gurobi_persistent import GurobiPersistent -from pyomo.contrib.solver.solvers.gurobi_direct import GurobiDirect +from pyomo.contrib.solver.solvers.gurobi.gurobi_direct import GurobiDirect +from pyomo.contrib.solver.solvers.gurobi.gurobi_persistent import ( + GurobiDirectQuadratic, + GurobiPersistent, +) from pyomo.contrib.solver.solvers.highs import Highs from pyomo.core.expr.numeric_expr import LinearExpression from pyomo.core.expr.compare import assertExpressionsEqual @@ -49,18 +52,27 @@ all_solvers = [ ('gurobi_persistent', GurobiPersistent), ('gurobi_direct', GurobiDirect), + ('gurobi_direct_quadratic', GurobiDirectQuadratic), ('ipopt', Ipopt), ('highs', Highs), ] mip_solvers = [ ('gurobi_persistent', GurobiPersistent), ('gurobi_direct', GurobiDirect), + ('gurobi_direct_quadratic', GurobiDirectQuadratic), ('highs', Highs), ] nlp_solvers = [('ipopt', Ipopt)] -qcp_solvers = [('gurobi_persistent', GurobiPersistent), ('ipopt', Ipopt)] +qcp_solvers = [ + ('gurobi_persistent', GurobiPersistent), + ('gurobi_direct_quadratic', GurobiDirectQuadratic), + ('ipopt', Ipopt), +] qp_solvers = qcp_solvers + [("highs", Highs)] -miqcqp_solvers = [('gurobi_persistent', GurobiPersistent)] +miqcqp_solvers = [ + ('gurobi_persistent', GurobiPersistent), + ('gurobi_direct_quadratic', GurobiDirectQuadratic), +] nl_solvers = [('ipopt', Ipopt)] nl_solvers_set = {i[0] for i in nl_solvers} @@ -1209,55 +1221,50 @@ def test_mutable_quadratic_objective_qp( def test_fixed_vars( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): - for treat_fixed_vars_as_params in [True, False]: - opt: SolverBase = opt_class() - if opt.is_persistent(): - opt = opt_class(treat_fixed_vars_as_params=treat_fixed_vars_as_params) - if not opt.available(): - raise unittest.SkipTest(f'Solver {opt.name} not available.') - if any(name.startswith(i) for i in nl_solvers_set): - if use_presolve: - opt.config.writer_config.linear_presolve = True - else: - opt.config.writer_config.linear_presolve = False - m = pyo.ConcreteModel() - m.x = pyo.Var() - m.x.fix(0) - m.y = pyo.Var() - a1 = 1 - a2 = -1 - b1 = 1 - b2 = 2 - m.obj = pyo.Objective(expr=m.y) - m.c1 = pyo.Constraint(expr=m.y >= a1 * m.x + b1) - m.c2 = pyo.Constraint(expr=m.y >= a2 * m.x + b2) - res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0) - self.assertAlmostEqual(m.y.value, 2) - m.x.unfix() - res = opt.solve(m) - self.assertAlmostEqual(m.x.value, (b2 - b1) / (a1 - a2)) - self.assertAlmostEqual(m.y.value, a1 * (b2 - b1) / (a1 - a2) + b1) - m.x.fix(0) - res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0) - self.assertAlmostEqual(m.y.value, 2) - m.x.value = 2 - res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 2) - self.assertAlmostEqual(m.y.value, 3) - m.x.value = 0 - res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0) - self.assertAlmostEqual(m.y.value, 2) + opt: SolverBase = opt_class() + if not opt.available(): + raise unittest.SkipTest(f'Solver {opt.name} not available.') + if any(name.startswith(i) for i in nl_solvers_set): + if use_presolve: + opt.config.writer_config.linear_presolve = True + else: + opt.config.writer_config.linear_presolve = False + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.x.fix(0) + m.y = pyo.Var() + a1 = 1 + a2 = -1 + b1 = 1 + b2 = 2 + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= a1 * m.x + b1) + m.c2 = pyo.Constraint(expr=m.y >= a2 * m.x + b2) + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 2) + m.x.unfix() + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, (b2 - b1) / (a1 - a2)) + self.assertAlmostEqual(m.y.value, a1 * (b2 - b1) / (a1 - a2) + b1) + m.x.fix(0) + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 2) + m.x.value = 2 + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 2) + self.assertAlmostEqual(m.y.value, 3) + m.x.value = 0 + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 2) @parameterized.expand(input=_load_tests(all_solvers)) def test_fixed_vars_2( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): opt: SolverBase = opt_class() - if opt.is_persistent(): - opt = opt_class(treat_fixed_vars_as_params=True) if not opt.available(): raise unittest.SkipTest(f'Solver {opt.name} not available.') if any(name.startswith(i) for i in nl_solvers_set): @@ -1301,8 +1308,6 @@ def test_fixed_vars_3( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): opt: SolverBase = opt_class() - if opt.is_persistent(): - opt = opt_class(treat_fixed_vars_as_params=True) if not opt.available(): raise unittest.SkipTest(f'Solver {opt.name} not available.') if any(name.startswith(i) for i in nl_solvers_set): @@ -1325,8 +1330,6 @@ def test_fixed_vars_4( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): opt: SolverBase = opt_class() - if opt.is_persistent(): - opt = opt_class(treat_fixed_vars_as_params=True) if not opt.available(): raise unittest.SkipTest(f'Solver {opt.name} not available.') if any(name.startswith(i) for i in nl_solvers_set): @@ -1880,10 +1883,7 @@ def test_fixed_binaries( res = opt.solve(m) self.assertAlmostEqual(res.incumbent_objective, 1) - if opt.is_persistent(): - opt: SolverBase = opt_class(treat_fixed_vars_as_params=False) - else: - opt = opt_class() + opt = opt_class() m.x.fix(0) res = opt.solve(m) self.assertAlmostEqual(res.incumbent_objective, 0) @@ -2037,33 +2037,30 @@ def test_bug_2(self, name: str, opt_class: Type[SolverBase], use_presolve: bool) This test is for a bug where an objective containing a fixed variable does not get updated properly when the variable is unfixed. """ - for fixed_var_option in [True, False]: - opt: SolverBase = opt_class() - if opt.is_persistent(): - opt = opt_class(treat_fixed_vars_as_params=fixed_var_option) - if not opt.available(): - raise unittest.SkipTest(f'Solver {opt.name} not available.') - if any(name.startswith(i) for i in nl_solvers_set): - if use_presolve: - opt.config.writer_config.linear_presolve = True - else: - opt.config.writer_config.linear_presolve = False - - m = pyo.ConcreteModel() - m.x = pyo.Var(bounds=(-10, 10)) - m.y = pyo.Var() - m.obj = pyo.Objective(expr=3 * m.y - m.x) - m.c = pyo.Constraint(expr=m.y >= m.x) - - m.x.fix(1) - res = opt.solve(m) - self.assertAlmostEqual(res.incumbent_objective, 2, 5) + opt: SolverBase = opt_class() + if not opt.available(): + raise unittest.SkipTest(f'Solver {opt.name} not available.') + if any(name.startswith(i) for i in nl_solvers_set): + if use_presolve: + opt.config.writer_config.linear_presolve = True + else: + opt.config.writer_config.linear_presolve = False - m.x.unfix() - m.x.setlb(-9) - m.x.setub(9) - res = opt.solve(m) - self.assertAlmostEqual(res.incumbent_objective, -18, 5) + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10)) + m.y = pyo.Var() + m.obj = pyo.Objective(expr=3 * m.y - m.x) + m.c = pyo.Constraint(expr=m.y >= m.x) + + m.x.fix(1) + res = opt.solve(m) + self.assertAlmostEqual(res.incumbent_objective, 2, 5) + + m.x.unfix() + m.x.setlb(-9) + m.x.setub(9) + res = opt.solve(m) + self.assertAlmostEqual(res.incumbent_objective, -18, 5) @parameterized.expand(input=_load_tests(nl_solvers)) def test_presolve_with_zero_coef(