# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright 2017 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.block import _BlockData, declare_custom_block import pyomo.environ as pyo from pyomo.solvers.plugins.solvers.persistent_solver import PersistentSolver from pyomo.core.expr.visitor import identify_variables from pyomo.common.collections import ComponentSet try: from mpi4py import MPI mpi4py_available = True except: mpi4py_available = False try: import numpy as np numpy_available = True except: numpy_available = False import logging logger = logging.getLogger(__name__) """ It is easier to understand this code after reading "A note on feasibility in Benders Decomposition" by Grothey et al. Original problem: min f(x, y) + h0(y) s.t. g(x, y) <= 0 h(y) <= 0 where y are the complicating variables. Reformulate to min h0(y) + eta s.t. g(x, y) <= 0 f(x, y) <= eta h(y) <= 0 Master problem must be of the form min h0(y) + eta s.t. h(y) <= 0 benders cuts where the last constraint will be generated automatically with BendersCutGenerators. The BendersCutGenerators must be handed a subproblem of the form min f(x, y) s.t. g(x, y) <= 0 except the constraints don't actually have to be in this form. The subproblem will automatically be transformed to min _z s.t. g(x, y) - z <= 0 (alpha) f(x, y) - eta - z <= 0 (beta) y - y_k = 0 (gamma) eta - eta_k = 0 (delta) """ solver_dual_sign_convention = dict() solver_dual_sign_convention['ipopt'] = -1 solver_dual_sign_convention['gurobi'] = -1 solver_dual_sign_convention['gurobi_direct'] = -1 solver_dual_sign_convention['gurobi_persistent'] = -1 solver_dual_sign_convention['cplex'] = -1 solver_dual_sign_convention['cplex_direct'] = -1 solver_dual_sign_convention['cplexdirect'] = -1 solver_dual_sign_convention['cplex_persistent'] = -1 solver_dual_sign_convention['glpk'] = -1 solver_dual_sign_convention['cbc'] = -1 solver_dual_sign_convention['xpress_direct'] = -1 solver_dual_sign_convention['xpress_persistent'] = -1 def _del_con(c): parent = c.parent_component() if parent.is_indexed(): parent.__delitem__(c.index()) else: assert parent is c c.parent_block().del_component(c) def _any_common_elements(a, b): if len(a) < len(b): for i in a: if i in b: return True else: for i in b: if i in a: return True return False def _setup_subproblem(b, master_vars, relax_subproblem_cons): # first get the objective and turn it into a constraint master_vars = ComponentSet(master_vars) objs = list(b.component_data_objects(pyo.Objective, descend_into=False, active=True)) if len(objs) != 1: raise ValueError('Subproblem must have exactly one objective') orig_obj = objs[0] orig_obj_expr = orig_obj.expr b.del_component(orig_obj) b._z = pyo.Var(bounds=(0, None)) b.objective = pyo.Objective(expr=b._z) b.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) b._eta = pyo.Var() b.aux_cons = pyo.ConstraintList() for c in list(b.component_data_objects(pyo.Constraint, descend_into=True, active=True, sort=True)): if not relax_subproblem_cons: c_vars = ComponentSet(identify_variables(c.body, include_fixed=False)) if not _any_common_elements(master_vars, c_vars): continue if c.equality: body = c.body rhs = pyo.value(c.lower) body -= rhs b.aux_cons.add(body - b._z <= 0) b.aux_cons.add(-body - b._z <= 0) _del_con(c) else: body = c.body lower = pyo.value(c.lower) upper = pyo.value(c.upper) if upper is not None: body_upper = body - upper - b._z b.aux_cons.add(body_upper <= 0) if lower is not None: body_lower = body - lower body_lower = -body_lower body_lower -= b._z b.aux_cons.add(body_lower <= 0) _del_con(c) b.obj_con = pyo.Constraint(expr=orig_obj_expr - b._eta - b._z <= 0) @declare_custom_block(name='BendersCutGenerator') class BendersCutGeneratorData(_BlockData): def __init__(self, component): if not mpi4py_available: raise ImportError('BendersCutGenerator requires mpi4py.') if not numpy_available: raise ImportError('BendersCutGenerator requires numpy.') _BlockData.__init__(self, component) self.num_subproblems_by_rank = 0 #np.zeros(self.comm.Get_size()) self.subproblems = list() self.complicating_vars_maps = list() self.master_vars = list() self.master_vars_indices = pyo.ComponentMap() self.master_etas = list() self.cuts = None self.subproblem_solvers = list() self.tol = None self.all_master_etas = list() self._subproblem_ndx_map = dict() # map from ndx in self.subproblems (local) to the global subproblem ndx def global_num_subproblems(self): return int(self.num_subproblems_by_rank.sum()) def local_num_subproblems(self): return len(self.subproblems) def set_input(self, master_vars, tol=1e-6, comm = None): """ It is very important for master_vars to be in the same order for every process. Parameters ---------- master_vars tol """ self.comm = None if comm is not None: self.comm = comm else: self.comm = MPI.COMM_WORLD self.num_subproblems_by_rank = np.zeros(self.comm.Get_size()) del self.cuts self.cuts = pyo.ConstraintList() self.subproblems = list() self.master_etas = list() self.complicating_vars_maps = list() self.master_vars = list(master_vars) self.master_vars_indices = pyo.ComponentMap() for i, v in enumerate(self.master_vars): self.master_vars_indices[v] = i self.tol = tol self.subproblem_solvers = list() self.all_master_etas = list() self._subproblem_ndx_map = dict() def add_subproblem(self, subproblem_fn, subproblem_fn_kwargs, master_eta, subproblem_solver='gurobi_persistent', relax_subproblem_cons=False): _rank = np.argmin(self.num_subproblems_by_rank) self.num_subproblems_by_rank[_rank] += 1 self.all_master_etas.append(master_eta) if _rank == self.comm.Get_rank(): self.master_etas.append(master_eta) subproblem, complicating_vars_map = subproblem_fn(**subproblem_fn_kwargs) self.subproblems.append(subproblem) self.complicating_vars_maps.append(complicating_vars_map) _setup_subproblem(subproblem, master_vars=[complicating_vars_map[i] for i in self.master_vars if i in complicating_vars_map], relax_subproblem_cons=relax_subproblem_cons) self._subproblem_ndx_map[len(self.subproblems) - 1] = self.global_num_subproblems() - 1 if isinstance(subproblem_solver, str): subproblem_solver = pyo.SolverFactory(subproblem_solver) self.subproblem_solvers.append(subproblem_solver) if isinstance(subproblem_solver, PersistentSolver): subproblem_solver.set_instance(subproblem) def generate_cut(self): coefficients = np.zeros(self.global_num_subproblems() * len(self.master_vars), dtype='d') constants = np.zeros(self.global_num_subproblems(), dtype='d') eta_coeffs = np.zeros(self.global_num_subproblems(), dtype='d') for local_subproblem_ndx in range(len(self.subproblems)): subproblem = self.subproblems[local_subproblem_ndx] global_subproblem_ndx = self._subproblem_ndx_map[local_subproblem_ndx] complicating_vars_map = self.complicating_vars_maps[local_subproblem_ndx] master_eta = self.master_etas[local_subproblem_ndx] coeff_ndx = global_subproblem_ndx * len(self.master_vars) subproblem.fix_complicating_vars = pyo.ConstraintList() var_to_con_map = pyo.ComponentMap() for master_var in self.master_vars: if master_var in complicating_vars_map: sub_var = complicating_vars_map[master_var] sub_var.value = master_var.value new_con = subproblem.fix_complicating_vars.add(sub_var - master_var.value == 0) var_to_con_map[master_var] = new_con subproblem.fix_eta = pyo.Constraint(expr=subproblem._eta - master_eta.value == 0) subproblem._eta.value = master_eta.value subproblem_solver = self.subproblem_solvers[local_subproblem_ndx] if subproblem_solver.name not in solver_dual_sign_convention: raise NotImplementedError('BendersCutGenerator is unaware of the dual sign convention of subproblem solver ' + subproblem_solver.name) sign_convention = solver_dual_sign_convention[subproblem_solver.name] if isinstance(subproblem_solver, PersistentSolver): for c in subproblem.fix_complicating_vars.values(): subproblem_solver.add_constraint(c) subproblem_solver.add_constraint(subproblem.fix_eta) res = subproblem_solver.solve(tee=False, load_solutions=False, save_results=False) if res.solver.termination_condition != pyo.TerminationCondition.optimal: raise RuntimeError('Unable to generate cut because subproblem failed to converge.') subproblem_solver.load_vars() subproblem_solver.load_duals() else: res = subproblem_solver.solve(subproblem, tee=False, load_solutions=False) if res.solver.termination_condition != pyo.TerminationCondition.optimal: raise RuntimeError('Unable to generate cut because subproblem failed to converge.') subproblem.solutions.load_from(res) constants[global_subproblem_ndx] = pyo.value(subproblem._z) eta_coeffs[global_subproblem_ndx] = sign_convention * pyo.value(subproblem.dual[subproblem.obj_con]) for master_var in self.master_vars: if master_var in complicating_vars_map: c = var_to_con_map[master_var] coefficients[coeff_ndx] = sign_convention * pyo.value(subproblem.dual[c]) coeff_ndx += 1 if isinstance(subproblem_solver, PersistentSolver): for c in subproblem.fix_complicating_vars.values(): subproblem_solver.remove_constraint(c) subproblem_solver.remove_constraint(subproblem.fix_eta) del subproblem.fix_complicating_vars del subproblem.fix_complicating_vars_index del subproblem.fix_eta total_num_subproblems = self.global_num_subproblems() global_constants = np.zeros(total_num_subproblems, dtype='d') global_coeffs = np.zeros(total_num_subproblems*len(self.master_vars), dtype='d') global_eta_coeffs = np.zeros(total_num_subproblems, dtype='d') comm = self.comm comm.Allreduce([constants, MPI.DOUBLE], [global_constants, MPI.DOUBLE]) comm.Allreduce([eta_coeffs, MPI.DOUBLE], [global_eta_coeffs, MPI.DOUBLE]) comm.Allreduce([coefficients, MPI.DOUBLE], [global_coeffs, MPI.DOUBLE]) global_constants = [float(i) for i in global_constants] global_coeffs = [float(i) for i in global_coeffs] global_eta_coeffs = [float(i) for i in global_eta_coeffs] coeff_ndx = 0 cuts_added = list() for global_subproblem_ndx in range(total_num_subproblems): cut_expr = global_constants[global_subproblem_ndx] if cut_expr > self.tol: master_eta = self.all_master_etas[global_subproblem_ndx] cut_expr -= global_eta_coeffs[global_subproblem_ndx] * (master_eta - master_eta.value) for master_var in self.master_vars: coeff = global_coeffs[coeff_ndx] cut_expr -= coeff * (master_var - master_var.value) coeff_ndx += 1 new_cut = self.cuts.add(cut_expr <= 0) cuts_added.append(new_cut) else: coeff_ndx += len(self.master_vars) return cuts_added