diff --git a/doc/OnlineDocs/developer_reference/solvers.rst b/doc/OnlineDocs/developer_reference/solvers.rst index 94fb684236f..9e3281246f4 100644 --- a/doc/OnlineDocs/developer_reference/solvers.rst +++ b/doc/OnlineDocs/developer_reference/solvers.rst @@ -45,9 +45,12 @@ with existing interfaces). * - Ipopt - ``ipopt`` - ``ipopt_v2`` - * - Gurobi + * - Gurobi (persistent) - ``gurobi`` - ``gurobi_v2`` + * - Gurobi (direct) + - ``gurobi_direct`` + - ``gurobi_direct_v2`` Using the new interfaces through the legacy interface ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/pyomo/contrib/solver/gurobi_direct.py b/pyomo/contrib/solver/gurobi_direct.py new file mode 100644 index 00000000000..edca7018f92 --- /dev/null +++ b/pyomo/contrib/solver/gurobi_direct.py @@ -0,0 +1,420 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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 datetime +import io +import math +import os + +from pyomo.common.config import ConfigValue +from pyomo.common.collections import ComponentMap, ComponentSet +from pyomo.common.dependencies import attempt_import +from pyomo.common.enums import ObjectiveSense +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.contrib.solver.base import SolverBase +from pyomo.contrib.solver.config import BranchAndBoundConfig +from pyomo.contrib.solver.results import Results, SolutionStatus, TerminationCondition +from pyomo.contrib.solver.solution import SolutionLoaderBase + +from pyomo.core.staleflag import StaleFlagManager + +from pyomo.repn.plugins.standard_form import LinearStandardFormCompiler + +gurobipy, gurobipy_available = attempt_import('gurobipy') + + +class GurobiConfig(BranchAndBoundConfig): + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + super(GurobiConfig, self).__init__( + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + 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._num_instances += 1 + + def __del__(self): + if not python_is_shutting_down(): + GurobiDirect._num_instances -= 1 + if GurobiDirect._num_instances == 0: + GurobiDirect.release_license() + + def load_vars(self, vars_to_load=None, solution_number=0): + assert solution_number == 0 + if self._grb_model.SolCount == 0: + raise RuntimeError( + 'Solver does not currently have a valid solution. Please ' + 'check the termination condition.' + ) + + 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 RuntimeError( + 'Solver does not currently have a valid solution. Please ' + 'check the termination condition.' + ) + + 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 RuntimeError( + 'Solver does not currently have valid duals. Please ' + 'check the termination condition.' + ) + + 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 RuntimeError( + 'Solver does not currently have valid reduced costs. Please ' + 'check the termination condition.' + ) + + 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(SolverBase): + CONFIG = GurobiConfig() + + _available = None + _num_instances = 0 + _tc_map = None + + def __init__(self, **kwds): + super().__init__(**kwds) + GurobiDirect._num_instances += 1 + + def available(self): + if not gurobipy_available: # this triggers the deferred import + return self.Availability.NotFound + elif self._available == self.Availability.BadVersion: + return self.Availability.BadVersion + else: + return self._check_license() + + def _check_license(self): + avail = False + try: + # Gurobipy writes out license file information when creating + # the environment + with capture_output(capture_fd=True): + m = gurobipy.Model() + avail = True + except gurobipy.GurobiError: + avail = False + + if avail: + if self._available is None: + self._available = GurobiDirect._check_full_license(m) + return self._available + else: + return self.Availability.BadLicense + + @classmethod + def _check_full_license(cls, model=None): + if model is None: + model = gurobipy.Model() + model.setParam('OutputFlag', 0) + try: + model.addVars(range(2001)) + model.optimize() + return cls.Availability.FullLicense + except gurobipy.GurobiError: + return cls.Availability.LimitedLicense + + def __del__(self): + if not python_is_shutting_down(): + GurobiDirect._num_instances -= 1 + if GurobiDirect._num_instances == 0: + self.release_license() + + @staticmethod + def release_license(): + if gurobipy_available: + with capture_output(capture_fd=True): + gurobipy.disposeDefaultEnv() + + def version(self): + version = ( + gurobipy.GRB.VERSION_MAJOR, + gurobipy.GRB.VERSION_MINOR, + gurobipy.GRB.VERSION_TECHNICAL, + ) + return version + + def solve(self, model, **kwds) -> Results: + start_timestamp = datetime.datetime.now(datetime.timezone.utc) + config = self.config(value=kwds, preserve_implicit=True) + 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') + + if len(repn.objectives) > 1: + raise ValueError( + 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 + lb = [] + ub = [] + for v in repn.columns: + _l, _u = v.bounds + if _l is None: + _l = ninf + if _u is None: + _u = inf + lb.append(_l) + ub.append(_u) + 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 = '=<>' # Note: ordering matches 0, 1, -1 + sense = [sense_type[r[1]] for r in repn.rows] + timer.stop('prepare_matrices') + + ostreams = [io.StringIO()] + config.tee + + try: + orig_cwd = os.getcwd() + if config.working_dir: + os.chdir(config.working_dir) + with TeeStream(*ostreams) as t, capture_output(t.STDOUT, capture_fd=False): + gurobi_model = gurobipy.Model() + + 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()) + timer.stop('transfer_model') + + options = config.solver_options + + gurobi_model.setParam('LogToConsole', 1) + + if config.threads is not None: + gurobi_model.setParam('Threads', config.threads) + if config.time_limit is not None: + gurobi_model.setParam('TimeLimit', config.time_limit) + if config.rel_gap is not None: + gurobi_model.setParam('MIPGap', config.rel_gap) + if config.abs_gap is not None: + gurobi_model.setParam('MIPGapAbs', config.abs_gap) + + if config.use_mipstart: + raise MouseTrap("MIPSTART not yet supported") + + for key, option in options.items(): + gurobi_model.setParam(key, option) + + timer.start('optimize') + gurobi_model.optimize() + timer.stop('optimize') + finally: + os.chdir(orig_cwd) + + res = self._postsolve( + timer, + config, + GurobiDirectSolutionLoader( + gurobi_model, A, x, repn.rows, repn.columns, repn.objectives + ), + ) + res.solver_configuration = 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 + status = grb_model.Status + + results = Results() + results.solution_loader = loader + results.timing_info.gurobi_time = grb_model.Runtime + + if grb_model.SolCount > 0: + if status == gurobipy.GRB.OPTIMAL: + results.solution_status = SolutionStatus.optimal + else: + results.solution_status = SolutionStatus.feasible + else: + results.solution_status = SolutionStatus.noSolution + + results.termination_condition = self._get_tc_map().get( + status, TerminationCondition.unknown + ) + + if ( + results.termination_condition + != TerminationCondition.convergenceCriteriaSatisfied + and config.raise_exception_on_nonoptimal_result + ): + raise RuntimeError( + 'Solver did not find the optimal solution. Set ' + 'opt.config.raise_exception_on_nonoptimal_result=False ' + 'to bypass this error.' + ) + + if loader._pyo_obj: + try: + if math.isfinite(grb_model.ObjVal): + results.incumbent_objective = grb_model.ObjVal + else: + results.incumbent_objective = None + except (gurobipy.GurobiError, AttributeError): + results.incumbent_objective = None + try: + results.objective_bound = grb_model.ObjBound + except (gurobipy.GurobiError, AttributeError): + if grb_model.ModelSense == ObjectiveSense.minimize: + results.objective_bound = -math.inf + else: + results.objective_bound = math.inf + else: + results.incumbent_objective = None + results.objective_bound = None + + results.iteration_count = grb_model.getAttr('IterCount') + + timer.start('load solution') + if config.load_solutions: + if grb_model.SolCount > 0: + results.solution_loader.load_vars() + else: + raise RuntimeError( + 'A feasible solution was not found, so no solution can be loaded.' + 'Please set opt.config.load_solutions=False and check ' + 'results.solution_status and ' + 'results.incumbent_objective before loading a solution.' + ) + timer.stop('load solution') + + return results + + 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 diff --git a/pyomo/contrib/solver/plugins.py b/pyomo/contrib/solver/plugins.py index 1a471d3bd06..82c10a32fd8 100644 --- a/pyomo/contrib/solver/plugins.py +++ b/pyomo/contrib/solver/plugins.py @@ -13,6 +13,7 @@ from .factory import SolverFactory from .ipopt import Ipopt from .gurobi import Gurobi +from .gurobi_direct import GurobiDirect def load(): @@ -22,3 +23,8 @@ def load(): SolverFactory.register( name='gurobi', legacy_name='gurobi_v2', doc='Persistent interface to Gurobi' )(Gurobi) + SolverFactory.register( + name='gurobi_direct', + legacy_name='gurobi_direct_v2', + doc='Direct (scipy-based) interface to Gurobi', + )(GurobiDirect) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index a4f4a3bc389..f91de2287b7 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -21,6 +21,7 @@ from pyomo.contrib.solver.base import SolverBase from pyomo.contrib.solver.ipopt import Ipopt from pyomo.contrib.solver.gurobi import Gurobi +from pyomo.contrib.solver.gurobi_direct import GurobiDirect from pyomo.core.expr.numeric_expr import LinearExpression @@ -32,8 +33,8 @@ if not param_available: raise unittest.SkipTest('Parameterized is not available.') -all_solvers = [('gurobi', Gurobi), ('ipopt', Ipopt)] -mip_solvers = [('gurobi', Gurobi)] +all_solvers = [('gurobi', Gurobi), ('gurobi_direct', GurobiDirect), ('ipopt', Ipopt)] +mip_solvers = [('gurobi', Gurobi), ('gurobi_direct', GurobiDirect)] nlp_solvers = [('ipopt', Ipopt)] qcp_solvers = [('gurobi', Gurobi), ('ipopt', Ipopt)] miqcqp_solvers = [('gurobi', Gurobi)] diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 315b2160d37..e684829e2f4 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -20,6 +20,7 @@ document_kwargs_from_configdict, ) from pyomo.common.dependencies import scipy, numpy as np +from pyomo.common.enums import ObjectiveSense from pyomo.common.gc_manager import PauseGC from pyomo.common.timing import TicTocTimer @@ -61,12 +62,16 @@ class LinearStandardFormInfo(object): Attributes ---------- - c : scipy.sparse.csr_array + c : scipy.sparse.csc_array The objective coefficients. Note that this is a sparse array and may contain multiple rows (for multiobjective problems). The objectives may be calculated by "c @ x" + c_offset : numpy.ndarray + + The list of objective constant offsets + A : scipy.sparse.csc_array The constraint coefficients. The constraint bodies may be @@ -89,6 +94,10 @@ class LinearStandardFormInfo(object): The list of Pyomo variable objects corresponding to columns in the `A` and `c` matrices. + objectives : List[ObjectiveData] + + The list of Pyomo objective objects corresponding to the active objectives + eliminated_vars: List[Tuple[VarData, NumericExpression]] The list of variables from the original model that do not appear @@ -101,12 +110,14 @@ class LinearStandardFormInfo(object): """ - def __init__(self, c, A, rhs, rows, columns, eliminated_vars): + def __init__(self, c, c_offset, A, rhs, rows, columns, objectives, eliminated_vars): self.c = c + self.c_offset = c_offset self.A = A self.rhs = rhs self.rows = rows self.columns = columns + self.objectives = objectives self.eliminated_vars = eliminated_vars @property @@ -148,6 +159,14 @@ class LinearStandardFormCompiler(object): 'mix of <=, ==, and >=)', ), ) + CONFIG.declare( + 'set_sense', + ConfigValue( + default=ObjectiveSense.minimize, + domain=InEnum(ObjectiveSense), + description='If not None, map all objectives to the specified sense.', + ), + ) CONFIG.declare( 'show_section_timing', ConfigValue( @@ -305,21 +324,19 @@ def write(self, model): # # Process objective # - if not component_map[Objective]: - objectives = [Objective(expr=1)] - objectives[0].construct() - else: - objectives = [] - for blk in component_map[Objective]: - objectives.extend( - blk.component_data_objects( - Objective, active=True, descend_into=False, sort=sorter - ) + set_sense = self.config.set_sense + objectives = [] + for blk in component_map[Objective]: + objectives.extend( + blk.component_data_objects( + Objective, active=True, descend_into=False, sort=sorter ) + ) + obj_offset = [] obj_data = [] obj_index = [] obj_index_ptr = [0] - for i, obj in enumerate(objectives): + for obj in objectives: repn = visitor.walk_expression(obj.expr) if repn.nonlinear is not None: raise ValueError( @@ -328,8 +345,10 @@ def write(self, model): ) N = len(repn.linear) obj_data.append(np.fromiter(repn.linear.values(), float, N)) - if obj.sense == maximize: + obj_offset.append(repn.constant) + if set_sense is not None and set_sense != obj.sense: obj_data[-1] *= -1 + obj_offset[-1] *= -1 obj_index.append( np.fromiter(map(var_order.__getitem__, repn.linear), float, N) ) @@ -456,13 +475,17 @@ def write(self, model): # Get the variable list columns = list(var_map.values()) # Convert the compiled data to scipy sparse matrices + if obj_data: + obj_data = np.concatenate(obj_data) + obj_index = np.concatenate(obj_index) c = scipy.sparse.csr_array( - (np.concatenate(obj_data), np.concatenate(obj_index), obj_index_ptr), - [len(obj_index_ptr) - 1, len(columns)], + (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, len(columns)] ).tocsc() + if rows: + con_data = np.concatenate(con_data) + con_index = np.concatenate(con_index) A = scipy.sparse.csr_array( - (np.concatenate(con_data), np.concatenate(con_index), con_index_ptr), - [len(rows), len(columns)], + (con_data, con_index, con_index_ptr), [len(rows), len(columns)] ).tocsc() # Some variables in the var_map may not actually appear in the @@ -495,7 +518,9 @@ def write(self, model): else: eliminated_vars = [] - info = LinearStandardFormInfo(c, A, rhs, rows, columns, eliminated_vars) + info = LinearStandardFormInfo( + c, np.array(obj_offset), A, rhs, rows, columns, objectives, eliminated_vars + ) timer.toc("Generated linear standard form representation", delta=False) return info