From 60376864685e0f9a09a3df3842a253bd2f80111a Mon Sep 17 00:00:00 2001 From: Qi Chen Date: Fri, 1 Feb 2019 11:36:51 -0500 Subject: [PATCH 1/6] Renaming GDPlbb to GDPbb --- .../contrib/{gdplbb/GDPlbb.py => gdpbb/GDPbb.py} | 16 ++++++++-------- pyomo/contrib/gdpbb/__init__.py | 1 + pyomo/contrib/gdpbb/plugins.py | 2 ++ .../test_GDPlbb.py => gdpbb/test_GDPbb.py} | 2 +- pyomo/contrib/gdplbb/__init__.py | 1 - pyomo/contrib/gdplbb/plugins.py | 2 -- pyomo/environ/__init__.py | 2 +- 7 files changed, 13 insertions(+), 13 deletions(-) rename pyomo/contrib/{gdplbb/GDPlbb.py => gdpbb/GDPbb.py} (96%) create mode 100644 pyomo/contrib/gdpbb/__init__.py create mode 100644 pyomo/contrib/gdpbb/plugins.py rename pyomo/contrib/{gdplbb/test_GDPlbb.py => gdpbb/test_GDPbb.py} (50%) delete mode 100644 pyomo/contrib/gdplbb/__init__.py delete mode 100644 pyomo/contrib/gdplbb/plugins.py diff --git a/pyomo/contrib/gdplbb/GDPlbb.py b/pyomo/contrib/gdpbb/GDPbb.py similarity index 96% rename from pyomo/contrib/gdplbb/GDPlbb.py rename to pyomo/contrib/gdpbb/GDPbb.py index 191857ca3ea..5b66d7777c9 100644 --- a/pyomo/contrib/gdplbb/GDPlbb.py +++ b/pyomo/contrib/gdpbb/GDPbb.py @@ -20,11 +20,11 @@ import heapq -@SolverFactory.register('gdplbb', +@SolverFactory.register('gdpbb', doc='Branch and Bound based GDP Solver') -class GDPlbbSolver(object): +class GDPbbSolver(object): """A branch and bound-based GDP solver.""" - CONFIG = ConfigBlock("gdplbb") + CONFIG = ConfigBlock("gdpbb") CONFIG.declare("solver", ConfigValue( default = "baron", description="Solver to use, defaults to baron" @@ -47,7 +47,7 @@ def solve(self, model, **kwds): config = self.CONFIG(kwds.pop('options',{})) config.set_value(kwds) - #Validate model to be used with gdplbb + #Validate model to be used with gdpbb self.validate_model(model) #Set solver as an MINLP solver = SolverFactory(config.solver) @@ -65,7 +65,7 @@ def solve(self, model, **kwds): objectives = model.component_data_objects(Objective, active=True) obj = next(objectives, None) obj_sign = 1 if obj.sense == minimize else -1 - print obj_sign + print(obj_sign) #clone original model for root node of branch and bound root = model.clone() @@ -105,7 +105,7 @@ def solve(self, model, **kwds): #initialize minheap for Branch and Bound algorithm heap = [] heapq.heappush(heap,(obj_sign * obj_value,root)) - print [i[0] for i in heap] + print( [i[0] for i in heap]) #loop to branch through the tree n = 0 while len(heap)>0: @@ -115,7 +115,7 @@ def solve(self, model, **kwds): #print [i[0] for i in heap] mdl = mdlpack[1] - print mdlpack[0] + print( mdlpack[0]) #if all the originally active disjunctions are active, solve and #return solution if(len(getattr(mdl,init_active_disjunctions_name)) == 0): @@ -182,7 +182,7 @@ def minlp_solve(self,gdp,solver,config): obj = next(objectives, None) obj_sign = 1 if obj.sense == minimize else -1 return obj_sign*float('inf') - delete(minlp) + def __enter__(self): return self diff --git a/pyomo/contrib/gdpbb/__init__.py b/pyomo/contrib/gdpbb/__init__.py new file mode 100644 index 00000000000..ddc637288f4 --- /dev/null +++ b/pyomo/contrib/gdpbb/__init__.py @@ -0,0 +1 @@ +import pyomo.contrib.gdpbb.GDPbb diff --git a/pyomo/contrib/gdpbb/plugins.py b/pyomo/contrib/gdpbb/plugins.py new file mode 100644 index 00000000000..5150afbfd2e --- /dev/null +++ b/pyomo/contrib/gdpbb/plugins.py @@ -0,0 +1,2 @@ +def load(): + import pyomo.contrib.gdpbb.GDPbb diff --git a/pyomo/contrib/gdplbb/test_GDPlbb.py b/pyomo/contrib/gdpbb/test_GDPbb.py similarity index 50% rename from pyomo/contrib/gdplbb/test_GDPlbb.py rename to pyomo/contrib/gdpbb/test_GDPbb.py index d3f5a12faa9..8b137891791 100644 --- a/pyomo/contrib/gdplbb/test_GDPlbb.py +++ b/pyomo/contrib/gdpbb/test_GDPbb.py @@ -1 +1 @@ - + diff --git a/pyomo/contrib/gdplbb/__init__.py b/pyomo/contrib/gdplbb/__init__.py deleted file mode 100644 index 8b35453fc1e..00000000000 --- a/pyomo/contrib/gdplbb/__init__.py +++ /dev/null @@ -1 +0,0 @@ -import pyomo.contrib.gdplbb.GDPlbb diff --git a/pyomo/contrib/gdplbb/plugins.py b/pyomo/contrib/gdplbb/plugins.py deleted file mode 100644 index 32f97bee962..00000000000 --- a/pyomo/contrib/gdplbb/plugins.py +++ /dev/null @@ -1,2 +0,0 @@ -def load(): - import pyomo.contrib.gdplbb.GDPlbb diff --git a/pyomo/environ/__init__.py b/pyomo/environ/__init__.py index 1c2da15faef..8efc4fc90e8 100644 --- a/pyomo/environ/__init__.py +++ b/pyomo/environ/__init__.py @@ -48,7 +48,7 @@ def _do_import(pkg_name): 'pyomo.contrib.example', 'pyomo.contrib.preprocessing', 'pyomo.contrib.gdpopt', - 'pyomo.contrib.gdplbb', + 'pyomo.contrib.gdpbb', 'pyomo.contrib.gdp_bounds', 'pyomo.contrib.trustregion', 'pyomo.contrib.multistart', From 8b227048a7b605ad2134a1d62eb6fed0be4cb250 Mon Sep 17 00:00:00 2001 From: Qi Chen Date: Fri, 1 Feb 2019 11:55:32 -0500 Subject: [PATCH 2/6] Add testing infrastructure --- pyomo/contrib/gdpbb/test_GDPbb.py | 74 +++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/pyomo/contrib/gdpbb/test_GDPbb.py b/pyomo/contrib/gdpbb/test_GDPbb.py index 8b137891791..c7f06bdaa60 100644 --- a/pyomo/contrib/gdpbb/test_GDPbb.py +++ b/pyomo/contrib/gdpbb/test_GDPbb.py @@ -1 +1,75 @@ +"""Tests for the GDPopt solver plugin.""" +import logging +from math import fabs +from os.path import abspath, dirname, join, normpath + +import pyutilib.th as unittest +from pyomo.environ import ConcreteModel, Objective, SolverFactory, Var, value, Integers, Block, Constraint, maximize +from pyomo.gdp import Disjunct, Disjunction +from pyutilib.misc import import_file +from pyomo.opt import TerminationCondition + +currdir = dirname(abspath(__file__)) +exdir = normpath(join(currdir, '..', '..', '..', 'examples', 'gdp')) + +minlp_solver = 'baron' + + +@unittest.skipIf(not SolverFactory(minlp_solver).available(), + "Required subsolver %s is not available" + % (minlp_solver,)) +class TestGLOA(unittest.TestCase): + """Tests for global logic-based outer approximation.""" + + def test_GLOA_8PP(self): + """Test the global logic-based outer approximation algorithm.""" + exfile = import_file( + join(exdir, 'eight_process', 'eight_proc_model.py')) + eight_process = exfile.build_eight_process_flowsheet() + SolverFactory('gdpbb').solve( + eight_process, tee=False, + solver=minlp_solver, + ) + self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2) + + def test_GLOA_strip_pack_default_init(self): + """Test logic-based outer approximation with strip packing.""" + exfile = import_file( + join(exdir, 'strip_packing', 'strip_packing_concrete.py')) + strip_pack = exfile.build_rect_strip_packing_model() + SolverFactory('gdpbb').solve( + strip_pack, tee=False, + solver=minlp_solver, + ) + self.assertTrue( + fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2) + + def test_GLOA_constrained_layout_default_init(self): + """Test LOA with constrained layout.""" + exfile = import_file( + join(exdir, 'constrained_layout', 'cons_layout_model.py')) + cons_layout = exfile.build_constrained_layout_model() + SolverFactory('gdpbb').solve( + cons_layout, tee=False, + solver=minlp_solver, + ) + objective_value = value(cons_layout.min_dist_cost.expr) + self.assertTrue( + fabs(objective_value - 41573) <= 200, + "Objective value of %s instead of 41573" % objective_value) + + def test_GLOA_ex_633_trespalacios(self): + """Test LOA with Francisco thesis example.""" + exfile = import_file(join(exdir, 'small_lit', 'ex_633_trespalacios.py')) + model = exfile.build_simple_nonconvex_gdp() + SolverFactory('gdpbb').solve( + model, tee=False, + solver=minlp_solver, + ) + objective_value = value(model.obj.expr) + self.assertAlmostEqual(objective_value, 4.46, 2) + + +if __name__ == '__main__': + unittest.main() From 8b8d3e85c51042fa736e9b0e6713043fc160b747 Mon Sep 17 00:00:00 2001 From: Qi Chen Date: Fri, 1 Feb 2019 11:56:04 -0500 Subject: [PATCH 3/6] IDE auto-format code. Adding tee option (unused for now) --- pyomo/contrib/gdpbb/GDPbb.py | 170 +++++++++++++++++------------------ 1 file changed, 81 insertions(+), 89 deletions(-) diff --git a/pyomo/contrib/gdpbb/GDPbb.py b/pyomo/contrib/gdpbb/GDPbb.py index 5b66d7777c9..c5646633e74 100644 --- a/pyomo/contrib/gdpbb/GDPbb.py +++ b/pyomo/contrib/gdpbb/GDPbb.py @@ -1,38 +1,33 @@ -import logging - -from six import iteritems +import heapq -import pyomo.util.plugin -from pyomo.core.base import expr as EXPR -from pyomo.core.base import (Block, Constraint, ConstraintList, Expression, - Objective, Set, Suffix, TransformationFactory, - Var, maximize, minimize, value) -from pyomo.core.base.block import generate_cuid_names -from pyomo.core.base.symbolic import differentiate +from pyomo.common.config import (ConfigBlock, ConfigValue) +from pyomo.common.modeling import unique_component_name +from pyomo.core.base import ( + Objective, TransformationFactory, + minimize, value) from pyomo.gdp import Disjunct, Disjunction +from pyomo.opt import SolverFactory, SolverStatus from pyomo.opt import TerminationCondition as tc -from pyomo.opt import SolutionStatus, SolverFactory, SolverStatus -from pyomo.opt.results import ProblemSense, SolverResults -from pyomo.common.config import (ConfigBlock, ConfigList, ConfigValue, In, - NonNegativeFloat, NonNegativeInt, - add_docstring_list) -from pyomo.common.modeling import unique_component_name -import heapq @SolverFactory.register('gdpbb', - doc='Branch and Bound based GDP Solver') + doc='Branch and Bound based GDP Solver') class GDPbbSolver(object): """A branch and bound-based GDP solver.""" CONFIG = ConfigBlock("gdpbb") CONFIG.declare("solver", ConfigValue( - default = "baron", + default="baron", description="Solver to use, defaults to baron" )) CONFIG.declare("solver_args", ConfigValue( default={}, description="Dictionary of keyword arguments to pass to the solver." )) + CONFIG.declare("tee", ConfigValue( + default=False, + domain=bool, + description="Flag to stream solver output to console." + )) def available(self, exception_flag=True): """Check if solver is available. @@ -44,117 +39,114 @@ def available(self, exception_flag=True): return True def solve(self, model, **kwds): - config = self.CONFIG(kwds.pop('options',{})) + config = self.CONFIG(kwds.pop('options', {})) config.set_value(kwds) - #Validate model to be used with gdpbb + # Validate model to be used with gdpbb self.validate_model(model) - #Set solver as an MINLP + # Set solver as an MINLP solver = SolverFactory(config.solver) - #Initialize ist containing indicator vars for reupdating model after solving - indicator_list_name = unique_component_name(model,"_indicator_list") + # Initialize ist containing indicator vars for reupdating model after solving + indicator_list_name = unique_component_name(model, "_indicator_list") indicator_vars = [] for disjunction in model.component_data_objects( - ctype = Disjunction, active=True): + ctype=Disjunction, active=True): for disjunct in disjunction.disjuncts: indicator_vars.append(disjunct.indicator_var) setattr(model, indicator_list_name, indicator_vars) - #get objective sense + # get objective sense objectives = model.component_data_objects(Objective, active=True) obj = next(objectives, None) obj_sign = 1 if obj.sense == minimize else -1 print(obj_sign) - #clone original model for root node of branch and bound + # clone original model for root node of branch and bound root = model.clone() - #set up lists to keep track of which disjunctions have been covered. + # set up lists to keep track of which disjunctions have been covered. - #this list keeps track of the original disjunctions that were active and are soon to be inactive - init_active_disjunctions_name = unique_component_name(root,"_init_active_disjunctions") + # this list keeps track of the original disjunctions that were active and are soon to be inactive + init_active_disjunctions_name = unique_component_name(root, "_init_active_disjunctions") init_active_disjunctions = list(root.component_data_objects( - ctype = Disjunction, active=True)) - setattr(root,init_active_disjunctions_name, init_active_disjunctions) + ctype=Disjunction, active=True)) + setattr(root, init_active_disjunctions_name, init_active_disjunctions) - #this list keeps track of the disjunctions that have been activated by the branch and bound - curr_active_disjunctions_name = unique_component_name(root,"_curr_active_disjunctions") + # this list keeps track of the disjunctions that have been activated by the branch and bound + curr_active_disjunctions_name = unique_component_name(root, "_curr_active_disjunctions") curr_active_disjunctions = [] - setattr(root,curr_active_disjunctions_name, curr_active_disjunctions) + setattr(root, curr_active_disjunctions_name, curr_active_disjunctions) - #deactivate all disjunctions in the model - #self.indicate(root) - for djn in getattr(root,init_active_disjunctions_name): + # deactivate all disjunctions in the model + # self.indicate(root) + for djn in getattr(root, init_active_disjunctions_name): djn.deactivate() - #Deactivate all disjuncts in model. To be reactivated with disjunction - #is reactivated. + # Deactivate all disjuncts in model. To be reactivated with disjunction + # is reactivated. for disj in root.component_data_objects(Disjunct, active=True): disj.deactivate() disj.indicator_var - for djn in getattr(root,init_active_disjunctions_name): + for djn in getattr(root, init_active_disjunctions_name): djn.disjuncts[0].indicator_var = 1 self.indicate(root) - #Satisfiability check would go here - - #solve the root node - obj_value = self.minlp_solve(root,solver,config) + # Satisfiability check would go here + # solve the root node + obj_value = self.minlp_solve(root, solver, config) - #initialize minheap for Branch and Bound algorithm + # initialize minheap for Branch and Bound algorithm heap = [] - heapq.heappush(heap,(obj_sign * obj_value,root)) - print( [i[0] for i in heap]) - #loop to branch through the tree + heapq.heappush(heap, (obj_sign * obj_value, root)) + print([i[0] for i in heap]) + # loop to branch through the tree n = 0 - while len(heap)>0: - n=n+1 - #pop best model off of heap + while len(heap) > 0: + n = n + 1 + # pop best model off of heap mdlpack = heapq.heappop(heap) - #print [i[0] for i in heap] + # print [i[0] for i in heap] mdl = mdlpack[1] - print( mdlpack[0]) - #if all the originally active disjunctions are active, solve and - #return solution - if(len(getattr(mdl,init_active_disjunctions_name)) == 0): + print(mdlpack[0]) + # if all the originally active disjunctions are active, solve and + # return solution + if (len(getattr(mdl, init_active_disjunctions_name)) == 0): orig_var_list = getattr(model, indicator_list_name) best_soln_var_list = getattr(mdl, indicator_list_name) - for orig_var, new_var in zip(orig_var_list,best_soln_var_list): + for orig_var, new_var in zip(orig_var_list, best_soln_var_list): if not orig_var.is_fixed(): orig_var.value = new_var.value TransformationFactory('gdp.fix_disjuncts').apply_to(model) - return solver.solve(model,**config.solver_args) + return solver.solve(model, **config.solver_args) - disjunction = getattr(mdl,init_active_disjunctions_name).pop(0) + disjunction = getattr(mdl, init_active_disjunctions_name).pop(0) for disj in list(disjunction.disjuncts): disj.activate() disjunction.activate() - getattr(mdl,curr_active_disjunctions_name).append(disjunction) + getattr(mdl, curr_active_disjunctions_name).append(disjunction) for disj in list(disjunction.disjuncts): disj.indicator_var = 0 for disj in list(disjunction.disjuncts): disj.indicator_var = 1 mnew = mdl.clone() disj.indicator_var = 0 - obj_value = self.minlp_solve(mnew,solver,config) - print[value(d.indicator_var) for d in mnew.component_data_objects(Disjunct, active=True)] - #self.indicate(mnew) - djn_left = len(getattr(mdl,init_active_disjunctions_name)) - ordering_tuple = (obj_sign * obj_value,djn_left,-n) - heapq.heappush(heap,(ordering_tuple,mnew)) - - - - def validate_model(self,model): - #Validates that model has only exclusive disjunctions + obj_value = self.minlp_solve(mnew, solver, config) + print([value(d.indicator_var) for d in mnew.component_data_objects(Disjunct, active=True)]) + # self.indicate(mnew) + djn_left = len(getattr(mdl, init_active_disjunctions_name)) + ordering_tuple = (obj_sign * obj_value, djn_left, -n) + heapq.heappush(heap, (ordering_tuple, mnew)) + + def validate_model(self, model): + # Validates that model has only exclusive disjunctions for d in model.component_data_objects( - ctype = Disjunction, active=True): - if(not d.xor): + ctype=Disjunction, active=True): + if (not d.xor): raise ValueError('GDPlbb unable to handle ' - 'non-exclusive disjunctions') + 'non-exclusive disjunctions') objectives = model.component_data_objects(Objective, active=True) obj = next(objectives, None) if next(objectives, None) is not None: @@ -164,24 +156,24 @@ def validate_model(self,model): raise RuntimeError( "GDP LBB solver is unable to handle model with no active objective.") - def minlp_solve(self,gdp,solver,config): + def minlp_solve(self, gdp, solver, config): minlp = gdp.clone() TransformationFactory('gdp.fix_disjuncts').apply_to(minlp) for disjunct in minlp.component_data_objects( - ctype = Disjunct, active=True): - disjunct.deactivate() #TODO this is HACK + ctype=Disjunct, active=True): + disjunct.deactivate() # TODO this is HACK - result = solver.solve(minlp,**config.solver_args) + result = solver.solve(minlp, **config.solver_args) if (result.solver.status is SolverStatus.ok and result.solver.termination_condition is tc.optimal): - objectives = minlp.component_data_objects(Objective, active=True) - obj = next(objectives, None) - return value(obj.expr) + objectives = minlp.component_data_objects(Objective, active=True) + obj = next(objectives, None) + return value(obj.expr) else: - objectives = minlp.component_data_objects(Objective, active=True) - obj = next(objectives, None) - obj_sign = 1 if obj.sense == minimize else -1 - return obj_sign*float('inf') + objectives = minlp.component_data_objects(Objective, active=True) + obj = next(objectives, None) + obj_sign = 1 if obj.sense == minimize else -1 + return obj_sign * float('inf') def __enter__(self): return self @@ -189,8 +181,8 @@ def __enter__(self): def __exit__(self, t, v, traceback): pass - def indicate(self,model): + def indicate(self, model): for disjunction in model.component_data_objects( - ctype = Disjunction): + ctype=Disjunction): for disjunct in disjunction.disjuncts: - print (disjunction.name,disjunct.name,value(disjunct.indicator_var)) + print(disjunction.name, disjunct.name, value(disjunct.indicator_var)) From 556b2e9aec92abc9a66a6b845529622bf54a265c Mon Sep 17 00:00:00 2001 From: Qi Chen Date: Sat, 2 Feb 2019 00:17:53 -0500 Subject: [PATCH 4/6] intermediate commit so that I can continue working on my laptop --- pyomo/contrib/gdpbb/GDPbb.py | 202 +++++++++++++++++++---------------- 1 file changed, 108 insertions(+), 94 deletions(-) diff --git a/pyomo/contrib/gdpbb/GDPbb.py b/pyomo/contrib/gdpbb/GDPbb.py index c5646633e74..07359e81294 100644 --- a/pyomo/contrib/gdpbb/GDPbb.py +++ b/pyomo/contrib/gdpbb/GDPbb.py @@ -2,6 +2,7 @@ from pyomo.common.config import (ConfigBlock, ConfigValue) from pyomo.common.modeling import unique_component_name +from pyomo.contrib.gdpopt.util import create_utility_block from pyomo.core.base import ( Objective, TransformationFactory, minimize, value) @@ -10,6 +11,10 @@ from pyomo.opt import TerminationCondition as tc +class GDPbbSolveData(object): + pass + + @SolverFactory.register('gdpbb', doc='Branch and Bound based GDP Solver') class GDPbbSolver(object): @@ -46,99 +51,108 @@ def solve(self, model, **kwds): self.validate_model(model) # Set solver as an MINLP solver = SolverFactory(config.solver) - - # Initialize ist containing indicator vars for reupdating model after solving - indicator_list_name = unique_component_name(model, "_indicator_list") - indicator_vars = [] - for disjunction in model.component_data_objects( - ctype=Disjunction, active=True): - for disjunct in disjunction.disjuncts: - indicator_vars.append(disjunct.indicator_var) - setattr(model, indicator_list_name, indicator_vars) - - # get objective sense - objectives = model.component_data_objects(Objective, active=True) - obj = next(objectives, None) - obj_sign = 1 if obj.sense == minimize else -1 - print(obj_sign) - - # clone original model for root node of branch and bound - root = model.clone() - - # set up lists to keep track of which disjunctions have been covered. - - # this list keeps track of the original disjunctions that were active and are soon to be inactive - init_active_disjunctions_name = unique_component_name(root, "_init_active_disjunctions") - init_active_disjunctions = list(root.component_data_objects( - ctype=Disjunction, active=True)) - setattr(root, init_active_disjunctions_name, init_active_disjunctions) - - # this list keeps track of the disjunctions that have been activated by the branch and bound - curr_active_disjunctions_name = unique_component_name(root, "_curr_active_disjunctions") - curr_active_disjunctions = [] - setattr(root, curr_active_disjunctions_name, curr_active_disjunctions) - - # deactivate all disjunctions in the model - # self.indicate(root) - for djn in getattr(root, init_active_disjunctions_name): - djn.deactivate() - # Deactivate all disjuncts in model. To be reactivated with disjunction - # is reactivated. - for disj in root.component_data_objects(Disjunct, active=True): - disj.deactivate() - disj.indicator_var - for djn in getattr(root, init_active_disjunctions_name): - djn.disjuncts[0].indicator_var = 1 - - self.indicate(root) - # Satisfiability check would go here - - # solve the root node - obj_value = self.minlp_solve(root, solver, config) - - # initialize minheap for Branch and Bound algorithm - heap = [] - heapq.heappush(heap, (obj_sign * obj_value, root)) - print([i[0] for i in heap]) - # loop to branch through the tree - n = 0 - while len(heap) > 0: - n = n + 1 - # pop best model off of heap - mdlpack = heapq.heappop(heap) - # print [i[0] for i in heap] - mdl = mdlpack[1] - - print(mdlpack[0]) - # if all the originally active disjunctions are active, solve and - # return solution - if (len(getattr(mdl, init_active_disjunctions_name)) == 0): - orig_var_list = getattr(model, indicator_list_name) - best_soln_var_list = getattr(mdl, indicator_list_name) - for orig_var, new_var in zip(orig_var_list, best_soln_var_list): - if not orig_var.is_fixed(): - orig_var.value = new_var.value - TransformationFactory('gdp.fix_disjuncts').apply_to(model) - - return solver.solve(model, **config.solver_args) - - disjunction = getattr(mdl, init_active_disjunctions_name).pop(0) - for disj in list(disjunction.disjuncts): - disj.activate() - disjunction.activate() - getattr(mdl, curr_active_disjunctions_name).append(disjunction) - for disj in list(disjunction.disjuncts): - disj.indicator_var = 0 - for disj in list(disjunction.disjuncts): - disj.indicator_var = 1 - mnew = mdl.clone() - disj.indicator_var = 0 - obj_value = self.minlp_solve(mnew, solver, config) - print([value(d.indicator_var) for d in mnew.component_data_objects(Disjunct, active=True)]) - # self.indicate(mnew) - djn_left = len(getattr(mdl, init_active_disjunctions_name)) - ordering_tuple = (obj_sign * obj_value, djn_left, -n) - heapq.heappush(heap, (ordering_tuple, mnew)) + solve_data = GDPbbSolveData() + + with create_utility_block(model, 'GDPbb_utils', solve_data): + # Initialize ist containing indicator vars for reupdating model after solving + indicator_list_name = unique_component_name(model, "_indicator_list") + indicator_vars = [] + for disjunction in model.component_data_objects( + ctype=Disjunction, active=True): + for disjunct in disjunction.disjuncts: + indicator_vars.append(disjunct.indicator_var) + setattr(model, indicator_list_name, indicator_vars) + + # get objective sense + objectives = model.component_data_objects(Objective, active=True) + obj = next(objectives, None) + obj_sign = 1 if obj.sense == minimize else -1 + print(obj_sign) + + # clone original model for root node of branch and bound + root = model.clone() + + # set up lists to keep track of which disjunctions have been covered. + + # this list keeps track of the original disjunctions that were active and are soon to be inactive + init_active_disjunctions_name = unique_component_name(root, "_init_active_disjunctions") + init_active_disjunctions = list(root.component_data_objects( + ctype=Disjunction, active=True)) + setattr(root, init_active_disjunctions_name, init_active_disjunctions) + + # this list keeps track of the disjunctions that have been activated by the branch and bound + curr_active_disjunctions_name = unique_component_name(root, "_curr_active_disjunctions") + curr_active_disjunctions = [] + setattr(root, curr_active_disjunctions_name, curr_active_disjunctions) + + # deactivate all disjunctions in the model + # self.indicate(root) + for djn in getattr(root, init_active_disjunctions_name): + djn.deactivate() + # Deactivate all disjuncts in model. To be reactivated with disjunction + # is reactivated. + for disj in root.component_data_objects(Disjunct, active=True): + disj._deactivate_without_fixing_indicator() + + # TODO [Qi]: what is this for? + for djn in getattr(root, init_active_disjunctions_name): + djn.disjuncts[0].indicator_var = 1 + + self.indicate(root) + # Satisfiability check would go here + + # solve the root node + obj_value = self.minlp_solve(root, solver, config) + + # initialize minheap for Branch and Bound algorithm + # Heap structure: (ordering tuple, model) + # Ordering tuple: (objective value, disjunctions_left, -counter) + # - select solutions with lower objective value, + # then fewer disjunctions left to explore (depth first), + # then more recently encountered (tiebreaker) + heap = [] + counter = 0 + disjunctions_left = len(getattr(root, init_active_disjunctions_name)) + heapq.heappush(heap, ((obj_sign * obj_value, disjunctions_left, counter), root)) + print([i[0] for i in heap]) + # loop to branch through the tree + while len(heap) > 0: + counter += 1 + # pop best model off of heap + mdlpack = heapq.heappop(heap) + # print [i[0] for i in heap] + mdl = mdlpack[1] + + print(mdlpack[0]) + # if all the originally active disjunctions are active, solve and + # return solution + if (len(getattr(mdl, init_active_disjunctions_name)) == 0): + orig_disjuncts = model.GDPbb_utils.disjunct_list + best_soln_disjuncts = mdl.GDPbb_utils.disjunct_list + for orig_disjunct, soln_disjunct in zip(orig_disjuncts, best_soln_disjuncts): + if not orig_disjunct.indicator_var.is_fixed(): + orig_disjunct.indicator_var.value = soln_disjunct.indicator_var.value + TransformationFactory('gdp.fix_disjuncts').apply_to(model) + + return solver.solve(model, **config.solver_args) + + disjunction = getattr(mdl, init_active_disjunctions_name).pop(0) + for disj in list(disjunction.disjuncts): + disj.activate() + disjunction.activate() + getattr(mdl, curr_active_disjunctions_name).append(disjunction) + for disj in list(disjunction.disjuncts): + disj.indicator_var = 0 + for disj in list(disjunction.disjuncts): + disj.indicator_var = 1 + mnew = mdl.clone() + disj.indicator_var = 0 + obj_value = self.minlp_solve(mnew, solver, config) + print([value(d.indicator_var) for d in mnew.component_data_objects(Disjunct, active=True)]) + # self.indicate(mnew) + djn_left = len(getattr(mdl, init_active_disjunctions_name)) + ordering_tuple = (obj_sign * obj_value, djn_left, -counter) + heapq.heappush(heap, (ordering_tuple, mnew)) def validate_model(self, model): # Validates that model has only exclusive disjunctions @@ -185,4 +199,4 @@ def indicate(self, model): for disjunction in model.component_data_objects( ctype=Disjunction): for disjunct in disjunction.disjuncts: - print(disjunction.name, disjunct.name, value(disjunct.indicator_var)) + print(disjunction.name, disjunct.name, disjunct.indicator_var.value) From dab9e1750ce18b819044d9e4aac4a728e20a7fb5 Mon Sep 17 00:00:00 2001 From: Qi Chen Date: Sat, 2 Feb 2019 01:42:37 -0500 Subject: [PATCH 5/6] Use GAMS BARON as the solver for now --- pyomo/contrib/gdpbb/test_GDPbb.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/gdpbb/test_GDPbb.py b/pyomo/contrib/gdpbb/test_GDPbb.py index c7f06bdaa60..d3d7cce6544 100644 --- a/pyomo/contrib/gdpbb/test_GDPbb.py +++ b/pyomo/contrib/gdpbb/test_GDPbb.py @@ -13,8 +13,8 @@ currdir = dirname(abspath(__file__)) exdir = normpath(join(currdir, '..', '..', '..', 'examples', 'gdp')) -minlp_solver = 'baron' - +minlp_solver = 'gams' +minlp_args=dict(solver='baron') @unittest.skipIf(not SolverFactory(minlp_solver).available(), "Required subsolver %s is not available" @@ -30,6 +30,7 @@ def test_GLOA_8PP(self): SolverFactory('gdpbb').solve( eight_process, tee=False, solver=minlp_solver, + solver_args=minlp_args, ) self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2) @@ -41,6 +42,7 @@ def test_GLOA_strip_pack_default_init(self): SolverFactory('gdpbb').solve( strip_pack, tee=False, solver=minlp_solver, + solver_args=minlp_args, ) self.assertTrue( fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2) @@ -53,6 +55,7 @@ def test_GLOA_constrained_layout_default_init(self): SolverFactory('gdpbb').solve( cons_layout, tee=False, solver=minlp_solver, + solver_args=minlp_args, ) objective_value = value(cons_layout.min_dist_cost.expr) self.assertTrue( @@ -66,6 +69,7 @@ def test_GLOA_ex_633_trespalacios(self): SolverFactory('gdpbb').solve( model, tee=False, solver=minlp_solver, + solver_args=minlp_args, ) objective_value = value(model.obj.expr) self.assertAlmostEqual(objective_value, 4.46, 2) From cc130651da9552b418ce0dd78f9a88db26af8341 Mon Sep 17 00:00:00 2001 From: Qi Chen Date: Sat, 2 Feb 2019 01:43:03 -0500 Subject: [PATCH 6/6] Clean up implementations --- pyomo/contrib/gdpbb/GDPbb.py | 97 +++++++++++++++++------------------- 1 file changed, 45 insertions(+), 52 deletions(-) diff --git a/pyomo/contrib/gdpbb/GDPbb.py b/pyomo/contrib/gdpbb/GDPbb.py index 07359e81294..7eb671489eb 100644 --- a/pyomo/contrib/gdpbb/GDPbb.py +++ b/pyomo/contrib/gdpbb/GDPbb.py @@ -67,7 +67,6 @@ def solve(self, model, **kwds): objectives = model.component_data_objects(Objective, active=True) obj = next(objectives, None) obj_sign = 1 if obj.sense == minimize else -1 - print(obj_sign) # clone original model for root node of branch and bound root = model.clone() @@ -75,34 +74,31 @@ def solve(self, model, **kwds): # set up lists to keep track of which disjunctions have been covered. # this list keeps track of the original disjunctions that were active and are soon to be inactive - init_active_disjunctions_name = unique_component_name(root, "_init_active_disjunctions") - init_active_disjunctions = list(root.component_data_objects( - ctype=Disjunction, active=True)) - setattr(root, init_active_disjunctions_name, init_active_disjunctions) + root.GDPbb_utils.unenforced_disjunctions = list( + disjunction for disjunction in root.GDPbb_utils.disjunction_list if disjunction.active + ) # this list keeps track of the disjunctions that have been activated by the branch and bound - curr_active_disjunctions_name = unique_component_name(root, "_curr_active_disjunctions") - curr_active_disjunctions = [] - setattr(root, curr_active_disjunctions_name, curr_active_disjunctions) + root.GDPbb_utils.curr_active_disjunctions = [] # deactivate all disjunctions in the model # self.indicate(root) - for djn in getattr(root, init_active_disjunctions_name): + for djn in root.GDPbb_utils.unenforced_disjunctions: djn.deactivate() - # Deactivate all disjuncts in model. To be reactivated with disjunction + # Deactivate all disjuncts in model. To be reactivated when disjunction # is reactivated. for disj in root.component_data_objects(Disjunct, active=True): disj._deactivate_without_fixing_indicator() - # TODO [Qi]: what is this for? - for djn in getattr(root, init_active_disjunctions_name): - djn.disjuncts[0].indicator_var = 1 + # # TODO [Qi]: what is this for? + # for djn in getattr(root, init_active_disjunctions_name): + # djn.disjuncts[0].indicator_var = 1 - self.indicate(root) + # self.indicate(root) # Satisfiability check would go here # solve the root node - obj_value = self.minlp_solve(root, solver, config) + obj_value, result = self.minlp_solve(root, solver, config) # initialize minheap for Branch and Bound algorithm # Heap structure: (ordering tuple, model) @@ -112,47 +108,45 @@ def solve(self, model, **kwds): # then more recently encountered (tiebreaker) heap = [] counter = 0 - disjunctions_left = len(getattr(root, init_active_disjunctions_name)) - heapq.heappush(heap, ((obj_sign * obj_value, disjunctions_left, counter), root)) + disjunctions_left = len(root.GDPbb_utils.unenforced_disjunctions) + heapq.heappush(heap, ((obj_sign * obj_value, disjunctions_left, -counter), root, result)) print([i[0] for i in heap]) # loop to branch through the tree while len(heap) > 0: - counter += 1 # pop best model off of heap - mdlpack = heapq.heappop(heap) + sort_tup, mdl, mdl_result = heapq.heappop(heap) + _, disjunctions_left, _ = sort_tup # print [i[0] for i in heap] - mdl = mdlpack[1] - print(mdlpack[0]) + print(sort_tup) # if all the originally active disjunctions are active, solve and # return solution - if (len(getattr(mdl, init_active_disjunctions_name)) == 0): - orig_disjuncts = model.GDPbb_utils.disjunct_list - best_soln_disjuncts = mdl.GDPbb_utils.disjunct_list - for orig_disjunct, soln_disjunct in zip(orig_disjuncts, best_soln_disjuncts): - if not orig_disjunct.indicator_var.is_fixed(): - orig_disjunct.indicator_var.value = soln_disjunct.indicator_var.value - TransformationFactory('gdp.fix_disjuncts').apply_to(model) - - return solver.solve(model, **config.solver_args) - - disjunction = getattr(mdl, init_active_disjunctions_name).pop(0) - for disj in list(disjunction.disjuncts): - disj.activate() - disjunction.activate() - getattr(mdl, curr_active_disjunctions_name).append(disjunction) - for disj in list(disjunction.disjuncts): - disj.indicator_var = 0 - for disj in list(disjunction.disjuncts): - disj.indicator_var = 1 + if disjunctions_left == 0: + # Model is solved. Copy over solution values. + for orig_var, soln_var in zip(model.GDPbb_utils.variable_list, mdl.GDPbb_utils.variable_list): + orig_var.value = soln_var.value + return mdl_result + + next_disjunction = mdl.GDPbb_utils.unenforced_disjunctions.pop(0) + next_disjunction.activate() + mdl.GDPbb_utils.curr_active_disjunctions.append(next_disjunction) + for disj in next_disjunction.disjuncts: + disj._activate_without_unfixing_indicator() + if not disj.indicator_var.fixed: + disj.indicator_var = 0 # initially set all indicator vars to zero + for disj in next_disjunction.disjuncts: + if not disj.indicator_var.fixed: + disj.indicator_var = 1 mnew = mdl.clone() - disj.indicator_var = 0 - obj_value = self.minlp_solve(mnew, solver, config) - print([value(d.indicator_var) for d in mnew.component_data_objects(Disjunct, active=True)]) + if not disj.indicator_var.fixed: + disj.indicator_var = 0 + obj_value, result = self.minlp_solve(mnew, solver, config) + # print([value(d.indicator_var) for d in mnew.component_data_objects(Disjunct, active=True)]) # self.indicate(mnew) - djn_left = len(getattr(mdl, init_active_disjunctions_name)) + counter += 1 + djn_left = len(mdl.GDPbb_utils.unenforced_disjunctions) ordering_tuple = (obj_sign * obj_value, djn_left, -counter) - heapq.heappush(heap, (ordering_tuple, mnew)) + heapq.heappush(heap, (ordering_tuple, root, result)) def validate_model(self, model): # Validates that model has only exclusive disjunctions @@ -178,16 +172,15 @@ def minlp_solve(self, gdp, solver, config): disjunct.deactivate() # TODO this is HACK result = solver.solve(minlp, **config.solver_args) + main_obj = next(minlp.component_data_objects(Objective, active=True)) + obj_sign = 1 if main_obj.sense == minimize else -1 if (result.solver.status is SolverStatus.ok and result.solver.termination_condition is tc.optimal): - objectives = minlp.component_data_objects(Objective, active=True) - obj = next(objectives, None) - return value(obj.expr) + return value(main_obj.expr), result + elif result.solver.termination_condition is tc.unbounded: + return obj_sign * float('-inf'), result else: - objectives = minlp.component_data_objects(Objective, active=True) - obj = next(objectives, None) - obj_sign = 1 if obj.sense == minimize else -1 - return obj_sign * float('inf') + return obj_sign * float('inf'), result def __enter__(self): return self