From 2e2cb24e74d8715df8d0bd605dc97984d40eba4f Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 6 Jun 2025 12:03:28 -0600 Subject: [PATCH 01/26] fix a bug in bigm_mixin --- pyomo/gdp/plugins/bigm_mixin.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/bigm_mixin.py b/pyomo/gdp/plugins/bigm_mixin.py index 4899551f918..a8cd9051ed4 100644 --- a/pyomo/gdp/plugins/bigm_mixin.py +++ b/pyomo/gdp/plugins/bigm_mixin.py @@ -9,12 +9,15 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ +import logging + from pyomo.gdp import GDP_Error from pyomo.common.collections import ComponentSet from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor import pyomo.contrib.fbbt.interval as interval from pyomo.core import Suffix +logger = logging.getLogger('pyomo.gdp.bigm') def _convert_M_to_tuple(M, constraint, disjunct=None): if not isinstance(M, (tuple, list)): @@ -23,7 +26,7 @@ def _convert_M_to_tuple(M, constraint, disjunct=None): else: try: M = (-M, M) - except: + except Exception: logger.error( "Error converting scalar M-value %s " "to (-M,M). Is %s not a numeric type?" % (M, type(M)) From 72003c4b4a17e37535034ba134d93a4948e794fd Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 6 Jun 2025 12:03:59 -0600 Subject: [PATCH 02/26] non-working attempt at threaded mbigm --- pyomo/gdp/plugins/multiple_bigm.py | 290 ++++++++++++++++++++++++++++- 1 file changed, 282 insertions(+), 8 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index b1e2581beca..19214c57dd1 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -57,6 +57,9 @@ from weakref import ref as weakref_ref +from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor, Future +import os + logger = logging.getLogger('pyomo.gdp.mbigm') _trusted_solvers = { @@ -71,6 +74,65 @@ 'highs', } +def _calc_M(constraint, other_disjunct, scratch_name, unsuccessful_message, upper): + + if not other_disjunct.active: + # This should not happen anymore because we don't deactivate infeasible disjuncts + # until we're done calculating all the Ms. It's a performance loss but what can + # you do? + print("This still happens???") + # return (0, False) + raise Exception("i think this should not happen anymore") + + # avoid mutating this while other threads are working on it + print(f"Before clone, other disjunct looks like this: {other_disjunct.pprint()}") + # other_disjunct = other_disjunct.clone() + scratch = getattr(other_disjunct, scratch_name) + print(f"{scratch.pprint()=}") + if upper: + scratch.obj.expr = constraint.body - constraint.upper + scratch.obj.sense = maximize + else: + scratch.obj.expr = constraint.body - constraint.lower + scratch.obj.sense = minimize + + print(f"Other disjunct looks like this: {other_disjunct.pprint()}") + + solver = SolverFactory('baron') + print("About to call solver") + results = solver.solve(other_disjunct, load_solutions=False) + if results.solver.termination_condition is TerminationCondition.infeasible: + # [2/18/24]: TODO: After the solver rewrite is complete, we will not + # need this check since we can actually determine from the + # termination condition whether or not the solver proved + # infeasibility or just terminated at local infeasiblity. For now, + # while this is not complete, it catches most of the solvers we + # trust, and, unless someone is so pathological as to *rename* an + # untrusted solver using a trusted solver name, it will never do the + # *wrong* thing. + if any(s in solver.name for s in _trusted_solvers): + logger.debug( + "Disjunct '%s' is infeasible, deactivating." % other_disjunct.name + ) + return (0, True) + else: + # This is a solver that might report + # 'infeasible' for local infeasibility, so we + # can't deactivate with confidence. To be + # conservative, we'll just complain about + # it. Post-solver-rewrite we will want to change + # this so that we check for 'proven_infeasible' + # and then we can abandon this hack + raise GDP_Error(unsuccessful_message) + elif results.solver.termination_condition is not TerminationCondition.optimal: + print("We had a solver error!") + raise GDP_Error(unsuccessful_message) + else: + print("We were optimal; loading sol") + other_disjunct.solutions.load_from(results) + M = value(scratch.obj.expr) + print(f"{M=}") + return (M, False) def Solver(val): if isinstance(val, str): @@ -217,6 +279,18 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): """, ), ) + CONFIG.declare( + 'threads', + ConfigValue( + default=None, + domain=int, + description="Number of worker threads to use when estimating M values", + doc=""" + If not specified, use up to the number of cores available, minus + one (but at least one). + """ + ) + ) transformation_name = 'mbigm' def __init__(self): @@ -271,19 +345,216 @@ def _apply_to_impl(self, instance, **kwds): gdp_tree = self._get_gdp_tree_from_targets(instance, targets) preprocessed_targets = gdp_tree.reverse_topological_sort() - for t in preprocessed_targets: - if t.ctype is Disjunction: - self._transform_disjunctionData( - t, - t.index(), - parent_disjunct=gdp_tree.parent(t), - root_disjunct=gdp_tree.root_disjunct(t), - ) + arg_Ms = self._config.bigM if self._config.bigM is not None else {} + self._transform_disjunctionDatas_threaded(preprocessed_targets, arg_Ms, gdp_tree) + + # for t in preprocessed_targets: + # if t.ctype is Disjunction: + # self._transform_disjunctionData( + # t, + # t.index(), + # parent_disjunct=gdp_tree.parent(t), + # root_disjunct=gdp_tree.root_disjunct(t), + # ) # issue warnings about anything that was in the bigM args dict that we # didn't use _warn_for_unused_bigM_args(self._config.bigM, self.used_args, logger) + def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp_tree): + # We want to do this one Disjunction at a time, but we also want to calculate the + # Ms in parallel. So we do a slightly convoluted dance of iterating the + # Disjunctions to get a list of M calculation jobs, then actually calculating Ms, + # then iterating Disjunctions again to actually transform the constraints. + + # To-do list in form (constraint, other_disj, unsucc_msg, upper: bool) + jobs = [] + # map Disjunction -> set of active Disjuncts + active_disjuncts = ComponentMap() + # set of Constraints processed during special handling of bound constraints; we + # will deactivate these, but not until we're done calculating Ms + transformed_constraints = set() + # Scratch Blocks to end up appearing on each Disjunct; they will be mutated + # several times before finished finding Ms + scratch_blocks = {} + # Ms ready for use + Ms = {} + # Disjunctions and their setup components. We will need to return to these later + disjunction_setup = {} + + for t in preprocessed_targets: + if t.ctype is Disjunction: + # contents of: _transform_disjunctionData() + # check for issues + if gdp_tree.root_disjunct(t) is not None: + # We do not support nested because, unlike in regular bigM, the + # constraints are not fully relaxed when the exactly-one constraint + # is not enforced. (For example, in this model: [1 <= x <= 3, [1 <= + # y <= 5] v [6 <= y <= 10]] v [5 <= x <= 10, 15 <= y <= 20]), we + # would need to put the relaxed inner-disjunction constraints on the + # parent Disjunct and process them again. This means the order in + # which we transformed Disjuncts would change the calculated M + # values. This is crazy, so we skip it. + raise GDP_Error( + "Found nested Disjunction '%s'. The multiple bigm " + "transformation does not support nested GDPs. " + "Please flatten the model before calling the " + "transformation" % t.name + ) + if not t.xor: + # This transformation assumes it can relax constraints assuming that + # another Disjunct is chosen. If it could be possible to choose both + # then that logic might fail. + raise GDP_Error( + "Cannot do multiple big-M reformulation for " + "Disjunction '%s' with OR constraint. " + "Must be an XOR!" % t.name + ) + + # start doing transformation + (transBlock, algebraic_constraint) = self._setup_transform_disjunctionData( + t, gdp_tree.root_disjunct(t) + ) + disjunction_setup[t] = (transBlock, algebraic_constraint) + + # already got arg_Ms + active_disjuncts[t] = ComponentSet(disj for disj in t.disjuncts if disj.active) + # this method returns the constraints transformed on this disjunct, + # use update because we are now saving these from all disjuncts for + # later + if self._config.reduce_bound_constraints: + transformed_constraints.update(self._transform_bound_constraints( + active_disjuncts[t], transBlock, arg_Ms + )) + # Get the jobs to calculate missing M values for this Disjunction. We can + # skip this if we are only doing bound constraints. + + if not self._config.only_mbigm_bound_constraints: + all_vars = list(self._get_all_var_objects(active_disjuncts[t])) + for disjunct, other_disjunct in itertools.product( + active_disjuncts[t], active_disjuncts[t] + ): + if disjunct is other_disjunct: + continue + elif id(other_disjunct) in scratch_blocks: + scratch = scratch_blocks[id(other_disjunct)] + else: + scratch = scratch_blocks[id(other_disjunct)] = Block() + other_disjunct.add_component( + unique_component_name(other_disjunct, "scratch"), scratch + ) + scratch.obj = Objective(expr=0) # placeholder, but I want to + # take the name before I add a + # bunch of random reference + # objects. + + # If the writers don't assume Vars are declared on the Block + # being solved, we won't need this! + for v in all_vars: + ref = Reference(v) + scratch.add_component(unique_component_name(scratch, v.name), ref) + scratch_name = scratch.local_name + for constraint in disjunct.component_data_objects( + Constraint, + active=True, + descend_into=Block, + sort=SortComponents.deterministic, + ): + if constraint in transformed_constraints: + continue + # First check args + if (constraint, other_disjunct) in arg_Ms: + (lower_M, upper_M) = _convert_M_to_tuple( + arg_Ms[constraint, other_disjunct], constraint, other_disjunct + ) + self.used_args[constraint, other_disjunct] = (lower_M, upper_M) + else: + (lower_M, upper_M) = (None, None) + unsuccessful_solve_msg = ( + "Unsuccessful solve to calculate M value to " + "relax constraint '%s' on Disjunct '%s' when " + "Disjunct '%s' is selected." + % (constraint.name, disjunct.name, other_disjunct.name) + ) + if constraint.lower is not None and lower_M is None: + jobs.append((constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, False)) + if constraint.upper is not None and upper_M is None: + jobs.append((constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, True)) + # Let's avoid touching arg_Ms and just have one source of + # truth for the current Ms + Ms[constraint, other_disjunct] = (lower_M, upper_M) + transBlock._mbm_values[constraint, other_disjunct] = (lower_M, upper_M) + # We are now outside of the DisjunctionDatas loop + # threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) + threads = 1 + running = {} + print(f"Executing {len(jobs)} jobs.") + with ThreadPoolExecutor(max_workers=threads, thread_name_prefix="mbigm_worker_thread_") as executor: + for constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, upper in jobs: + running[constraint, other_disjunct, upper] = ( + executor.submit( + _calc_M, + constraint, + other_disjunct, + scratch_name, + unsuccessful_solve_msg, + upper, + ) + ) + # # cheat for testing purposes + # f = Future() + # f.set_result(_calc_M(constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, upper)) + # running[constraint, other_disjunct, upper] = f + + to_deactivate = [] + for (constraint, other_disjunct, upper), future in running.items(): + M, disjunct_infeasible = future.result() + print(f"Got result: {M}, {disjunct_infeasible}") + if disjunct_infeasible: + # If we made it here without an exception, the solver is on the + # trusted solvers list + logger.debug( + "Disjunct '%s' is infeasible; will deactivate." % other_disjunct.name + ) + # Don't deactivate live because we have concurrent jobs still running + if other_disjunct not in to_deactivate: + to_deactivate.append(other_disjunct) + continue + # Note that we can't just transform immediately because we might be + # waiting on the other one of upper_M/lower_M and I'd like to do + # constraints all at once. + if upper: + Ms[constraint, other_disjunct] = (Ms[constraint, other_disjunct][0], M) + else: + Ms[constraint, other_disjunct] = (M, Ms[constraint, other_disjunct][1]) + transBlock._mbm_values[constraint, other_disjunct] = Ms[constraint, other_disjunct] + + # No longer in threaded context + for disjunct in to_deactivate: + print("Deactivating a disjunct.") + disjunct.deactivate() + active_disjuncts[disjunct.parent_block()].remove(disjunct) + # We now have all the M values. Do some cleanup: + for con in transformed_constraints: + con.deactivate() + for block in scratch_blocks.values(): + block.parent_block().del_component(block) + + print(f"Final {Ms.values()=}") + # Iterate the Disjunctions again to actually transform them + for disjunction, (transBlock, algebraic_constraint) in disjunction_setup.items(): + or_expr = 0 + for disjunct in active_disjuncts[disjunction]: + or_expr += disjunct.indicator_var.get_associated_binary() + self._transform_disjunct(disjunct, transBlock, active_disjuncts[disjunction], Ms) + algebraic_constraint.add(disjunction.index(), or_expr == 1) + # map the DisjunctionData to its XOR constraint to mark it as + # transformed + disjunction._algebraic_constraint = weakref_ref(algebraic_constraint[disjunction.index()]) + disjunction.deactivate() + + + def _transform_disjunctionData(self, obj, index, parent_disjunct, root_disjunct): if root_disjunct is not None: # We do not support nested because, unlike in regular bigM, the @@ -681,6 +952,7 @@ def _solve_disjunct_for_M( solver = self._config.solver + print("About to call a solver") results = solver.solve(other_disjunct, load_solutions=False) if results.solver.termination_condition is TerminationCondition.infeasible: # [2/18/24]: TODO: After the solver rewrite is complete, we will not @@ -710,8 +982,10 @@ def _solve_disjunct_for_M( elif results.solver.termination_condition is not TerminationCondition.optimal: raise GDP_Error(unsuccessful_solve_msg) else: + print("We were optimal; loading sol") other_disjunct.solutions.load_from(results) M = value(scratch_block.obj.expr) + print(f"{M=}") return M def _warn_for_active_suffix(self, suffix, disjunct, active_disjuncts, Ms): From 526bca4e857b3a2b16cad98678094197cf0eb701 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 9 Jun 2025 12:27:59 -0600 Subject: [PATCH 03/26] Working attempt at parallel mbigm when using baron, but leads to a deadlock with gurobi --- pyomo/gdp/plugins/multiple_bigm.py | 322 +++++++++++------------------ 1 file changed, 125 insertions(+), 197 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 19214c57dd1..d4a87c0f1a5 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -37,10 +37,12 @@ Suffix, value, Var, + ConcreteModel, ) from pyomo.core.base import Reference, TransformationFactory import pyomo.core.expr as EXPR from pyomo.core.util import target_list +from pyomo.common.errors import DeveloperError from pyomo.gdp import Disjunct, Disjunction, GDP_Error from pyomo.gdp.plugins.bigm_mixin import ( @@ -74,65 +76,6 @@ 'highs', } -def _calc_M(constraint, other_disjunct, scratch_name, unsuccessful_message, upper): - - if not other_disjunct.active: - # This should not happen anymore because we don't deactivate infeasible disjuncts - # until we're done calculating all the Ms. It's a performance loss but what can - # you do? - print("This still happens???") - # return (0, False) - raise Exception("i think this should not happen anymore") - - # avoid mutating this while other threads are working on it - print(f"Before clone, other disjunct looks like this: {other_disjunct.pprint()}") - # other_disjunct = other_disjunct.clone() - scratch = getattr(other_disjunct, scratch_name) - print(f"{scratch.pprint()=}") - if upper: - scratch.obj.expr = constraint.body - constraint.upper - scratch.obj.sense = maximize - else: - scratch.obj.expr = constraint.body - constraint.lower - scratch.obj.sense = minimize - - print(f"Other disjunct looks like this: {other_disjunct.pprint()}") - - solver = SolverFactory('baron') - print("About to call solver") - results = solver.solve(other_disjunct, load_solutions=False) - if results.solver.termination_condition is TerminationCondition.infeasible: - # [2/18/24]: TODO: After the solver rewrite is complete, we will not - # need this check since we can actually determine from the - # termination condition whether or not the solver proved - # infeasibility or just terminated at local infeasiblity. For now, - # while this is not complete, it catches most of the solvers we - # trust, and, unless someone is so pathological as to *rename* an - # untrusted solver using a trusted solver name, it will never do the - # *wrong* thing. - if any(s in solver.name for s in _trusted_solvers): - logger.debug( - "Disjunct '%s' is infeasible, deactivating." % other_disjunct.name - ) - return (0, True) - else: - # This is a solver that might report - # 'infeasible' for local infeasibility, so we - # can't deactivate with confidence. To be - # conservative, we'll just complain about - # it. Post-solver-rewrite we will want to change - # this so that we check for 'proven_infeasible' - # and then we can abandon this hack - raise GDP_Error(unsuccessful_message) - elif results.solver.termination_condition is not TerminationCondition.optimal: - print("We had a solver error!") - raise GDP_Error(unsuccessful_message) - else: - print("We were optimal; loading sol") - other_disjunct.solutions.load_from(results) - M = value(scratch.obj.expr) - print(f"{M=}") - return (M, False) def Solver(val): if isinstance(val, str): @@ -348,35 +291,23 @@ def _apply_to_impl(self, instance, **kwds): arg_Ms = self._config.bigM if self._config.bigM is not None else {} self._transform_disjunctionDatas_threaded(preprocessed_targets, arg_Ms, gdp_tree) - # for t in preprocessed_targets: - # if t.ctype is Disjunction: - # self._transform_disjunctionData( - # t, - # t.index(), - # parent_disjunct=gdp_tree.parent(t), - # root_disjunct=gdp_tree.root_disjunct(t), - # ) - # issue warnings about anything that was in the bigM args dict that we # didn't use _warn_for_unused_bigM_args(self._config.bigM, self.used_args, logger) def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp_tree): - # We want to do this one Disjunction at a time, but we also want to calculate the - # Ms in parallel. So we do a slightly convoluted dance of iterating the - # Disjunctions to get a list of M calculation jobs, then actually calculating Ms, - # then iterating Disjunctions again to actually transform the constraints. + # We wish we could do this one Disjunction at a time, but we also want to + # calculate the Ms in parallel. So we do a slightly convoluted dance of iterating + # the Disjunctions to get a list of M calculation jobs, then actually calculating + # Ms, then iterating Disjunctions again to actually transform the constraints. - # To-do list in form (constraint, other_disj, unsucc_msg, upper: bool) + # To-do list in form (constraint, other_disj, scratch, unsucc_msg, upper) jobs = [] - # map Disjunction -> set of active Disjuncts + # map Disjunction -> set of its active Disjuncts active_disjuncts = ComponentMap() # set of Constraints processed during special handling of bound constraints; we # will deactivate these, but not until we're done calculating Ms transformed_constraints = set() - # Scratch Blocks to end up appearing on each Disjunct; they will be mutated - # several times before finished finding Ms - scratch_blocks = {} # Ms ready for use Ms = {} # Disjunctions and their setup components. We will need to return to these later @@ -384,8 +315,6 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp for t in preprocessed_targets: if t.ctype is Disjunction: - # contents of: _transform_disjunctionData() - # check for issues if gdp_tree.root_disjunct(t) is not None: # We do not support nested because, unlike in regular bigM, the # constraints are not fully relaxed when the exactly-one constraint @@ -436,24 +365,7 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp ): if disjunct is other_disjunct: continue - elif id(other_disjunct) in scratch_blocks: - scratch = scratch_blocks[id(other_disjunct)] - else: - scratch = scratch_blocks[id(other_disjunct)] = Block() - other_disjunct.add_component( - unique_component_name(other_disjunct, "scratch"), scratch - ) - scratch.obj = Objective(expr=0) # placeholder, but I want to - # take the name before I add a - # bunch of random reference - # objects. - - # If the writers don't assume Vars are declared on the Block - # being solved, we won't need this! - for v in all_vars: - ref = Reference(v) - scratch.add_component(unique_component_name(scratch, v.name), ref) - scratch_name = scratch.local_name + for constraint in disjunct.component_data_objects( Constraint, active=True, @@ -464,11 +376,13 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp continue # First check args if (constraint, other_disjunct) in arg_Ms: + print("Getting Ms from arg") (lower_M, upper_M) = _convert_M_to_tuple( arg_Ms[constraint, other_disjunct], constraint, other_disjunct ) self.used_args[constraint, other_disjunct] = (lower_M, upper_M) else: + print("Did not find Ms in arg") (lower_M, upper_M) = (None, None) unsuccessful_solve_msg = ( "Unsuccessful solve to calculate M value to " @@ -476,49 +390,66 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp "Disjunct '%s' is selected." % (constraint.name, disjunct.name, other_disjunct.name) ) + if constraint.lower is not None and lower_M is None: - jobs.append((constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, False)) + scratch = self._get_scratch_block(constraint, other_disjunct, all_vars, False) + jobs.append((constraint, other_disjunct, scratch, unsuccessful_solve_msg, False)) if constraint.upper is not None and upper_M is None: - jobs.append((constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, True)) + scratch = self._get_scratch_block(constraint, other_disjunct, all_vars, True) + jobs.append((constraint, other_disjunct, scratch, unsuccessful_solve_msg, True)) # Let's avoid touching arg_Ms and just have one source of # truth for the current Ms Ms[constraint, other_disjunct] = (lower_M, upper_M) transBlock._mbm_values[constraint, other_disjunct] = (lower_M, upper_M) - # We are now outside of the DisjunctionDatas loop - # threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) - threads = 1 + # (Now exiting: if not only_mbigm_bound_constraints) + # + # In this case we just pass through arg_Ms as-is and let + # transform_constraint parse it + # + # TODO revisit what's going on here. Can we just skip half the + # transformation instead of doing this? We have queued zero jobs. + else: + Ms = arg_Ms + # (Now exiting the DisjunctionDatas loop) + threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) + # threads = self._config.threads if self._config.threads is not None else 1 running = {} - print(f"Executing {len(jobs)} jobs.") - with ThreadPoolExecutor(max_workers=threads, thread_name_prefix="mbigm_worker_thread_") as executor: - for constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, upper in jobs: + print(f"Before processing, {Ms=}") + print(f"Executing {len(jobs)} jobs on {threads} worker threads.") + with ThreadPoolExecutor(max_workers=threads, thread_name_prefix="mBigM-Worker-Thread") as executor: + for constraint, other_disjunct, scratch, unsuccessful_solve_msg, upper in jobs: running[constraint, other_disjunct, upper] = ( executor.submit( - _calc_M, - constraint, - other_disjunct, - scratch_name, + self._calc_M, + scratch, unsuccessful_solve_msg, upper, ) ) # # cheat for testing purposes # f = Future() - # f.set_result(_calc_M(constraint, other_disjunct, scratch_name, unsuccessful_solve_msg, upper)) + # f.set_result(self._calc_M(scratch, unsuccessful_solve_msg, upper)) # running[constraint, other_disjunct, upper] = f - to_deactivate = [] + to_deactivate = set() for (constraint, other_disjunct, upper), future in running.items(): - M, disjunct_infeasible = future.result() - print(f"Got result: {M}, {disjunct_infeasible}") + try: + M, disjunct_infeasible = future.result() + print(f"Got result: {M}, {disjunct_infeasible}") + except Exception as e: + import traceback + print("@@@@@@@@@@@ PRINTING ERROR! @@@@@@@@@") + traceback.print_tb(e.__traceback__) + raise e if disjunct_infeasible: - # If we made it here without an exception, the solver is on the - # trusted solvers list - logger.debug( - "Disjunct '%s' is infeasible; will deactivate." % other_disjunct.name - ) # Don't deactivate live because we have concurrent jobs still running if other_disjunct not in to_deactivate: - to_deactivate.append(other_disjunct) + # If we made it here without an exception, the solver is on the + # trusted solvers list + logger.debug( + "Disjunct '%s' is infeasible, deactivating." % other_disjunct.name + ) + to_deactivate.add(other_disjunct) continue # Note that we can't just transform immediately because we might be # waiting on the other one of upper_M/lower_M and I'd like to do @@ -529,16 +460,24 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp Ms[constraint, other_disjunct] = (M, Ms[constraint, other_disjunct][1]) transBlock._mbm_values[constraint, other_disjunct] = Ms[constraint, other_disjunct] - # No longer in threaded context - for disjunct in to_deactivate: - print("Deactivating a disjunct.") - disjunct.deactivate() - active_disjuncts[disjunct.parent_block()].remove(disjunct) + # Now exiting threaded context + # TODO: make this not awful + if to_deactivate: + for disjunction in disjunction_setup.keys(): + for disjunct in set(active_disjuncts[disjunction]): + if disjunct in to_deactivate: + disjunct.deactivate() + active_disjuncts[disjunction].remove(disjunct) + # for disjunct in to_deactivate: + # print("Deactivating a disjunct.") + # disjunct.deactivate() + # want something like + # active_disjuncts[disjunct.parent_disjunction()].remove(disjunct), but there + # is no such method on Disjunct. Instead we will do something terrible + # We now have all the M values. Do some cleanup: for con in transformed_constraints: con.deactivate() - for block in scratch_blocks.values(): - block.parent_block().del_component(block) print(f"Final {Ms.values()=}") # Iterate the Disjunctions again to actually transform them @@ -553,78 +492,67 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp disjunction._algebraic_constraint = weakref_ref(algebraic_constraint[disjunction.index()]) disjunction.deactivate() + def _get_scratch_block(self, constraint, other_disjunct, all_vars, upper): + # Should this be a Block or a ConcreteModel? + scratch = Block() + if upper: + scratch.obj = Objective(expr=constraint.body - constraint.upper, sense=maximize) + else: + scratch.obj = Objective(expr=constraint.body - constraint.lower, sense=minimize) + # Create a Block component (reference) that actually points to a DisjunctData, as + # we want the writer to write it as if it were an ordinary Block rather than + # getting any special handling. DisjunctData inherits BlockData, so this should + # be legal. + scratch.disjunct_ref = Reference(other_disjunct, ctype=Block) + # If the writers don't assume Vars are declared on the Block + # being solved, we won't need this! + for v in all_vars: + ref = Reference(v) + scratch.add_component(unique_component_name(scratch, v.name), ref) + scratch.construct() # should this go in _calc_M? + return scratch + + def _calc_M(self, scratch, unsuccessful_message, upper): + + # scratch.pprint() + solver = SolverFactory('baron') + # solver = self._config.solver + print("About to call solver") + results = solver.solve(scratch, tee=False, load_solutions=False) - - def _transform_disjunctionData(self, obj, index, parent_disjunct, root_disjunct): - if root_disjunct is not None: - # We do not support nested because, unlike in regular bigM, the - # constraints are not fully relaxed when the exactly-one constraint - # is not enforced. (For example, in this model: [1 <= x <= 3, [1 <= - # y <= 5] v [6 <= y <= 10]] v [5 <= x <= 10, 15 <= y <= 20]), we - # would need to put the relaxed inner-disjunction constraints on the - # parent Disjunct and process them again. This means the order in - # which we transformed Disjuncts would change the calculated M - # values. This is crazy, so we skip it. - raise GDP_Error( - "Found nested Disjunction '%s'. The multiple bigm " - "transformation does not support nested GDPs. " - "Please flatten the model before calling the " - "transformation" % obj.name - ) - - if not obj.xor: - # This transformation assumes it can relax constraints assuming that - # another Disjunct is chosen. If it could be possible to choose both - # then that logic might fail. - raise GDP_Error( - "Cannot do multiple big-M reformulation for " - "Disjunction '%s' with OR constraint. " - "Must be an XOR!" % obj.name - ) - - (transBlock, algebraic_constraint) = self._setup_transform_disjunctionData( - obj, root_disjunct - ) - - ## Here's the logic for the actual transformation - - arg_Ms = self._config.bigM if self._config.bigM is not None else {} - - # ESJ: I am relying on the fact that the ComponentSet is going to be - # ordered here, but using a set because I will remove infeasible - # Disjuncts from it if I encounter them calculating M's. - active_disjuncts = ComponentSet(disj for disj in obj.disjuncts if disj.active) - # First handle the bound constraints if we are dealing with them - # separately - transformed_constraints = set() - if self._config.reduce_bound_constraints: - transformed_constraints = self._transform_bound_constraints( - active_disjuncts, transBlock, arg_Ms - ) - - Ms = arg_Ms - if not self._config.only_mbigm_bound_constraints: - Ms = transBlock.calculated_missing_m_values = ( - self._calculate_missing_M_values( - active_disjuncts, arg_Ms, transBlock, transformed_constraints - ) - ) - - # Now we can deactivate the constraints we deferred, so that we don't - # re-transform them - for cons in transformed_constraints: - cons.deactivate() - - or_expr = 0 - for disjunct in active_disjuncts: - or_expr += disjunct.indicator_var.get_associated_binary() - self._transform_disjunct(disjunct, transBlock, active_disjuncts, Ms) - algebraic_constraint.add(index, or_expr == 1) - # map the DisjunctionData to its XOR constraint to mark it as - # transformed - obj._algebraic_constraint = weakref_ref(algebraic_constraint[index]) - - obj.deactivate() + if results.solver.termination_condition is TerminationCondition.infeasible: + # [2/18/24]: TODO: After the solver rewrite is complete, we will not + # need this check since we can actually determine from the + # termination condition whether or not the solver proved + # infeasibility or just terminated at local infeasiblity. For now, + # while this is not complete, it catches most of the solvers we + # trust, and, unless someone is so pathological as to *rename* an + # untrusted solver using a trusted solver name, it will never do the + # *wrong* thing. + if any(s in solver.name for s in _trusted_solvers): + return (0, True) + else: + # This is a solver that might report + # 'infeasible' for local infeasibility, so we + # can't deactivate with confidence. To be + # conservative, we'll just complain about + # it. Post-solver-rewrite we will want to change + # this so that we check for 'proven_infeasible' + # and then we can abandon this hack + raise GDP_Error(unsuccessful_message) + elif results.solver.termination_condition is not TerminationCondition.optimal: + print("We had a solver error!") + raise GDP_Error(unsuccessful_message) + else: + print("We were optimal; getting M from objective bound") + # note: This tranformation can be made faster by allowing the solver a gap. As + # long as we have a bound, it's still valid (though not as tight) + if upper: + M = results.problem.upper_bound + else: + M = results.problem.lower_bound + print(f"{M=}") + return (M, False) def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms): # We've already filtered out deactivated disjuncts, so we know obj is From 69b22ceaa3eebbdf9806cb8716b692003ca606ab Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 9 Jun 2025 14:51:42 -0600 Subject: [PATCH 04/26] delete unused code --- pyomo/gdp/plugins/multiple_bigm.py | 148 ++--------------------------- 1 file changed, 9 insertions(+), 139 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index d4a87c0f1a5..f04558b5b8c 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -376,13 +376,13 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp continue # First check args if (constraint, other_disjunct) in arg_Ms: - print("Getting Ms from arg") + # print("Getting Ms from arg") (lower_M, upper_M) = _convert_M_to_tuple( arg_Ms[constraint, other_disjunct], constraint, other_disjunct ) self.used_args[constraint, other_disjunct] = (lower_M, upper_M) else: - print("Did not find Ms in arg") + # print("Did not find Ms in arg") (lower_M, upper_M) = (None, None) unsuccessful_solve_msg = ( "Unsuccessful solve to calculate M value to " @@ -414,7 +414,7 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) # threads = self._config.threads if self._config.threads is not None else 1 running = {} - print(f"Before processing, {Ms=}") + # print(f"Before processing, {Ms=}") print(f"Executing {len(jobs)} jobs on {threads} worker threads.") with ThreadPoolExecutor(max_workers=threads, thread_name_prefix="mBigM-Worker-Thread") as executor: for constraint, other_disjunct, scratch, unsuccessful_solve_msg, upper in jobs: @@ -435,7 +435,7 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp for (constraint, other_disjunct, upper), future in running.items(): try: M, disjunct_infeasible = future.result() - print(f"Got result: {M}, {disjunct_infeasible}") + # print(f"Got result: {M}, {disjunct_infeasible}") except Exception as e: import traceback print("@@@@@@@@@@@ PRINTING ERROR! @@@@@@@@@") @@ -479,7 +479,7 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp for con in transformed_constraints: con.deactivate() - print(f"Final {Ms.values()=}") + # print(f"Final {Ms.values()=}") # Iterate the Disjunctions again to actually transform them for disjunction, (transBlock, algebraic_constraint) in disjunction_setup.items(): or_expr = 0 @@ -517,7 +517,7 @@ def _calc_M(self, scratch, unsuccessful_message, upper): # scratch.pprint() solver = SolverFactory('baron') # solver = self._config.solver - print("About to call solver") + # print("About to call solver") results = solver.solve(scratch, tee=False, load_solutions=False) if results.solver.termination_condition is TerminationCondition.infeasible: @@ -541,17 +541,17 @@ def _calc_M(self, scratch, unsuccessful_message, upper): # and then we can abandon this hack raise GDP_Error(unsuccessful_message) elif results.solver.termination_condition is not TerminationCondition.optimal: - print("We had a solver error!") + # print("We had a solver error!") raise GDP_Error(unsuccessful_message) else: - print("We were optimal; getting M from objective bound") + # print("We were optimal; getting M from objective bound") # note: This tranformation can be made faster by allowing the solver a gap. As # long as we have a bound, it's still valid (though not as tight) if upper: M = results.problem.upper_bound else: M = results.problem.lower_bound - print(f"{M=}") + # print(f"{M=}") return (M, False) def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms): @@ -786,136 +786,6 @@ def _get_all_var_objects(self, active_disjuncts): seen.add(id(var)) yield var - def _calculate_missing_M_values( - self, active_disjuncts, arg_Ms, transBlock, transformed_constraints - ): - scratch_blocks = {} - all_vars = list(self._get_all_var_objects(active_disjuncts)) - for disjunct, other_disjunct in itertools.product( - active_disjuncts, active_disjuncts - ): - if disjunct is other_disjunct: - continue - elif id(other_disjunct) in scratch_blocks: - scratch = scratch_blocks[id(other_disjunct)] - else: - scratch = scratch_blocks[id(other_disjunct)] = Block() - other_disjunct.add_component( - unique_component_name(other_disjunct, "scratch"), scratch - ) - scratch.obj = Objective(expr=0) # placeholder, but I want to - # take the name before I add a - # bunch of random reference - # objects. - - # If the writers don't assume Vars are declared on the Block - # being solved, we won't need this! - for v in all_vars: - ref = Reference(v) - scratch.add_component(unique_component_name(scratch, v.name), ref) - - for constraint in disjunct.component_data_objects( - Constraint, - active=True, - descend_into=Block, - sort=SortComponents.deterministic, - ): - if constraint in transformed_constraints: - continue - # First check args - if (constraint, other_disjunct) in arg_Ms: - (lower_M, upper_M) = _convert_M_to_tuple( - arg_Ms[constraint, other_disjunct], constraint, other_disjunct - ) - self.used_args[constraint, other_disjunct] = (lower_M, upper_M) - else: - (lower_M, upper_M) = (None, None) - unsuccessful_solve_msg = ( - "Unsuccessful solve to calculate M value to " - "relax constraint '%s' on Disjunct '%s' when " - "Disjunct '%s' is selected." - % (constraint.name, disjunct.name, other_disjunct.name) - ) - if constraint.lower is not None and lower_M is None: - # last resort: calculate - if lower_M is None: - scratch.obj.expr = constraint.body - constraint.lower - scratch.obj.sense = minimize - lower_M = self._solve_disjunct_for_M( - other_disjunct, - scratch, - unsuccessful_solve_msg, - active_disjuncts, - ) - if constraint.upper is not None and upper_M is None: - # last resort: calculate - if upper_M is None: - scratch.obj.expr = constraint.body - constraint.upper - scratch.obj.sense = maximize - upper_M = self._solve_disjunct_for_M( - other_disjunct, - scratch, - unsuccessful_solve_msg, - active_disjuncts, - ) - arg_Ms[constraint, other_disjunct] = (lower_M, upper_M) - transBlock._mbm_values[constraint, other_disjunct] = (lower_M, upper_M) - - # clean up the scratch blocks - for blk in scratch_blocks.values(): - blk.parent_block().del_component(blk) - - return arg_Ms - - def _solve_disjunct_for_M( - self, other_disjunct, scratch_block, unsuccessful_solve_msg, active_disjuncts - ): - if not other_disjunct.active: - # If a Disjunct is infeasible, we will discover that and deactivate - # it when we are calculating the M values. We remove that disjunct - # from active_disjuncts inside of the loop in - # _calculate_missing_M_values. So that means that we might have - # deactivated Disjuncts here that we should skip over. - return 0 - - solver = self._config.solver - - print("About to call a solver") - results = solver.solve(other_disjunct, load_solutions=False) - if results.solver.termination_condition is TerminationCondition.infeasible: - # [2/18/24]: TODO: After the solver rewrite is complete, we will not - # need this check since we can actually determine from the - # termination condition whether or not the solver proved - # infeasibility or just terminated at local infeasiblity. For now, - # while this is not complete, it catches most of the solvers we - # trust, and, unless someone is so pathological as to *rename* an - # untrusted solver using a trusted solver name, it will never do the - # *wrong* thing. - if any(s in solver.name for s in _trusted_solvers): - logger.debug( - "Disjunct '%s' is infeasible, deactivating." % other_disjunct.name - ) - other_disjunct.deactivate() - active_disjuncts.remove(other_disjunct) - M = 0 - else: - # This is a solver that might report - # 'infeasible' for local infeasibility, so we - # can't deactivate with confidence. To be - # conservative, we'll just complain about - # it. Post-solver-rewrite we will want to change - # this so that we check for 'proven_infeasible' - # and then we can abandon this hack - raise GDP_Error(unsuccessful_solve_msg) - elif results.solver.termination_condition is not TerminationCondition.optimal: - raise GDP_Error(unsuccessful_solve_msg) - else: - print("We were optimal; loading sol") - other_disjunct.solutions.load_from(results) - M = value(scratch_block.obj.expr) - print(f"{M=}") - return M - def _warn_for_active_suffix(self, suffix, disjunct, active_disjuncts, Ms): if suffix.local_name == 'BigM': logger.debug( From cf09d528a251ee2f47023cd0242ad3e855bda5d6 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 10 Jun 2025 16:35:28 -0600 Subject: [PATCH 05/26] Improve several edge cases --- pyomo/gdp/plugins/multiple_bigm.py | 110 ++++++++++++++++++----------- pyomo/gdp/tests/test_mbigm.py | 2 +- 2 files changed, 68 insertions(+), 44 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index f04558b5b8c..434de7e945c 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -11,6 +11,7 @@ import itertools import logging +import math from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.common.config import ConfigDict, ConfigValue @@ -59,6 +60,7 @@ from weakref import ref as weakref_ref +import threading from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor, Future import os @@ -76,6 +78,8 @@ 'highs', } +# Keep Solver objects in thread-local storage +_thread_local = threading.local() def Solver(val): if isinstance(val, str): @@ -140,7 +144,7 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): CONFIG.declare( 'solver', ConfigValue( - default='gurobi', + default='baron', # TODO restore this to gurobi after fixing deadlock domain=Solver, description="A solver to use to solve the continuous subproblems for " "calculating the M values", @@ -229,11 +233,25 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): domain=int, description="Number of worker threads to use when estimating M values", doc=""" - If not specified, use up to the number of cores available, minus + If not specified, use up to the number of available cores, minus one (but at least one). """ ) ) + CONFIG.declare( + 'use_primal_bound', + ConfigValue( + default=False, + domain=bool, + description="When estimating M values, use the primal bound instead of dual bound", + doc=""" + This is necessary when using a local solver such as ipopt, but be + aware that interior feasible points for this subproblem do not give + valid values for M. That is, in the presence of any numerical error, + this option will lead to a slightly wrong reformulation. + """ + ) + ) transformation_name = 'mbigm' def __init__(self): @@ -298,19 +316,21 @@ def _apply_to_impl(self, instance, **kwds): def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp_tree): # We wish we could do this one Disjunction at a time, but we also want to # calculate the Ms in parallel. So we do a slightly convoluted dance of iterating - # the Disjunctions to get a list of M calculation jobs, then actually calculating - # Ms, then iterating Disjunctions again to actually transform the constraints. + # the Disjunctions to get a list of M calculation jobs, then calculating Ms + # multithreaded, then iterating Disjunctions again to actually transform the + # constraints. # To-do list in form (constraint, other_disj, scratch, unsucc_msg, upper) jobs = [] # map Disjunction -> set of its active Disjuncts active_disjuncts = ComponentMap() - # set of Constraints processed during special handling of bound constraints; we + # set of Constraints processed during special handling of bound constraints: we # will deactivate these, but not until we're done calculating Ms transformed_constraints = set() - # Ms ready for use + # Finished M values Ms = {} - # Disjunctions and their setup components. We will need to return to these later + # Disjunctions and their setup components. We will return to these after + # calculating Ms. disjunction_setup = {} for t in preprocessed_targets: @@ -357,7 +377,6 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp )) # Get the jobs to calculate missing M values for this Disjunction. We can # skip this if we are only doing bound constraints. - if not self._config.only_mbigm_bound_constraints: all_vars = list(self._get_all_var_objects(active_disjuncts[t])) for disjunct, other_disjunct in itertools.product( @@ -407,12 +426,12 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp # transform_constraint parse it # # TODO revisit what's going on here. Can we just skip half the - # transformation instead of doing this? We have queued zero jobs. + # transformation instead of doing this? We have queued zero jobs so + # none of the threads stuff will run. else: Ms = arg_Ms # (Now exiting the DisjunctionDatas loop) threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) - # threads = self._config.threads if self._config.threads is not None else 1 running = {} # print(f"Before processing, {Ms=}") print(f"Executing {len(jobs)} jobs on {threads} worker threads.") @@ -433,14 +452,14 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp to_deactivate = set() for (constraint, other_disjunct, upper), future in running.items(): - try: - M, disjunct_infeasible = future.result() - # print(f"Got result: {M}, {disjunct_infeasible}") - except Exception as e: - import traceback - print("@@@@@@@@@@@ PRINTING ERROR! @@@@@@@@@") - traceback.print_tb(e.__traceback__) - raise e + # try: + M, disjunct_infeasible = future.result() + # print(f"Got result: {M}, {disjunct_infeasible}") + # except Exception as e: + # import traceback + # print("@@@@@@@@@@@ PRINTING ERROR! @@@@@@@@@") + # traceback.print_tb(e.__traceback__) + # raise e if disjunct_infeasible: # Don't deactivate live because we have concurrent jobs still running if other_disjunct not in to_deactivate: @@ -452,30 +471,18 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp to_deactivate.add(other_disjunct) continue # Note that we can't just transform immediately because we might be - # waiting on the other one of upper_M/lower_M and I'd like to do - # constraints all at once. + # waiting on the other one of upper_M or lower_M. if upper: Ms[constraint, other_disjunct] = (Ms[constraint, other_disjunct][0], M) else: Ms[constraint, other_disjunct] = (M, Ms[constraint, other_disjunct][1]) transBlock._mbm_values[constraint, other_disjunct] = Ms[constraint, other_disjunct] - # Now exiting threaded context - # TODO: make this not awful - if to_deactivate: - for disjunction in disjunction_setup.keys(): - for disjunct in set(active_disjuncts[disjunction]): - if disjunct in to_deactivate: - disjunct.deactivate() - active_disjuncts[disjunction].remove(disjunct) - # for disjunct in to_deactivate: - # print("Deactivating a disjunct.") - # disjunct.deactivate() - # want something like - # active_disjuncts[disjunct.parent_disjunction()].remove(disjunct), but there - # is no such method on Disjunct. Instead we will do something terrible - - # We now have all the M values. Do some cleanup: + # Now exiting threaded context. Cleanup: + for disjunct in to_deactivate: + disjunct.deactivate() + active_disjuncts[gdp_tree.parent(disjunct)].remove(disjunct) + for con in transformed_constraints: con.deactivate() @@ -502,7 +509,7 @@ def _get_scratch_block(self, constraint, other_disjunct, all_vars, upper): # Create a Block component (reference) that actually points to a DisjunctData, as # we want the writer to write it as if it were an ordinary Block rather than # getting any special handling. DisjunctData inherits BlockData, so this should - # be legal. + # be valid. scratch.disjunct_ref = Reference(other_disjunct, ctype=Block) # If the writers don't assume Vars are declared on the Block # being solved, we won't need this! @@ -513,10 +520,13 @@ def _get_scratch_block(self, constraint, other_disjunct, all_vars, upper): return scratch def _calc_M(self, scratch, unsuccessful_message, upper): - - # scratch.pprint() - solver = SolverFactory('baron') - # solver = self._config.solver + # Retrieve a solver object from thread-local storage; if we don't have one, make + # a new solver reasonably similar to the one passed via config. + if hasattr(_thread_local, 'solver'): + solver = _thread_local.solver + else: + solver_arg = self._config.solver + solver = _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) # print("About to call solver") results = solver.solve(scratch, tee=False, load_solutions=False) @@ -546,11 +556,25 @@ def _calc_M(self, scratch, unsuccessful_message, upper): else: # print("We were optimal; getting M from objective bound") # note: This tranformation can be made faster by allowing the solver a gap. As - # long as we have a bound, it's still valid (though not as tight) - if upper: + # long as we have a bound, it's still valid (though not as tight). + + # We should use the dual bound here. The primal bound is mathematically + # incorrect in the presence of error, but local solvers need it. + if upper != self._config.use_primal_bound: M = results.problem.upper_bound else: M = results.problem.lower_bound + if not math.isfinite(M): + raise GDP_Error( + "Subproblem solved to optimality, but no finite dual bound was " + "obtained. If you would like to instead use the best obtained " + "solution, set the option use_primal_bound=True. This is " + "necessary when using a local solver such as ipopt, but be " + "aware that interior feasible points for this subproblem do " + "not give valid values for M." + if not self._config.use_primal_bound else + "No finite objective after optimal solve." + ) # print(f"{M=}") return (M, False) diff --git a/pyomo/gdp/tests/test_mbigm.py b/pyomo/gdp/tests/test_mbigm.py index aaa4750681c..56daf21e77f 100644 --- a/pyomo/gdp/tests/test_mbigm.py +++ b/pyomo/gdp/tests/test_mbigm.py @@ -1062,7 +1062,7 @@ def test_calculate_Ms_infeasible_Disjunct_local_solver(self): r"Disjunct 'disjunction_disjuncts\[0\]' is selected.", ): TransformationFactory('gdp.mbigm').apply_to( - m, solver=SolverFactory('ipopt'), reduce_bound_constraints=False + m, solver=SolverFactory('ipopt'), reduce_bound_constraints=False, use_primal_bound=True ) @unittest.skipUnless(gurobi_available, "Gurobi is not available") From 89218cf79b2d2e2ee72817bf4564715198fe1979 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 11 Jun 2025 16:07:05 -0600 Subject: [PATCH 06/26] Handle edge cases "better" --- pyomo/gdp/plugins/multiple_bigm.py | 95 ++++++++++++++++++++++++------ 1 file changed, 76 insertions(+), 19 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 434de7e945c..948f6b9d081 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -78,6 +78,16 @@ 'highs', } +# Solvers that do not report objective bounds or values in the result and therefore +# require loading the solution. For now, this prevents the use of multithreading; see the +# comment in _calc_M for more information. + +# TODO: this no longer prevents multithreading unless I go back to the singlethreaded +# method +_solvers_require_loading_solutions = { + 'ipopt', +} + # Keep Solver objects in thread-local storage _thread_local = threading.local() @@ -320,7 +330,7 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp # multithreaded, then iterating Disjunctions again to actually transform the # constraints. - # To-do list in form (constraint, other_disj, scratch, unsucc_msg, upper) + # To-do list in form (constraint, other_disj, scratch, unsucc_msg, is_upper) jobs = [] # map Disjunction -> set of its active Disjuncts active_disjuncts = ComponentMap() @@ -432,26 +442,31 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp Ms = arg_Ms # (Now exiting the DisjunctionDatas loop) threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) + ## Not doing this, unless I decide to do it again + # solver_name = self._config.solver.name + # if threads > 1 and solver_name in _solvers_require_loading_solutions: + # logger.warning(f"Solver '{solver_name}' cannot be used in multithreaded M calculations for now - reverting to a single thread. Pass the option threads=1 to silence this warning. ") + # threads = 1 running = {} - # print(f"Before processing, {Ms=}") - print(f"Executing {len(jobs)} jobs on {threads} worker threads.") + if jobs: + logger.info(f"Executing {len(jobs)} jobs on {threads} worker threads.") with ThreadPoolExecutor(max_workers=threads, thread_name_prefix="mBigM-Worker-Thread") as executor: - for constraint, other_disjunct, scratch, unsuccessful_solve_msg, upper in jobs: - running[constraint, other_disjunct, upper] = ( + for constraint, other_disjunct, scratch, unsuccessful_solve_msg, is_upper in jobs: + running[constraint, other_disjunct, is_upper] = ( executor.submit( self._calc_M, scratch, unsuccessful_solve_msg, - upper, + is_upper, ) ) # # cheat for testing purposes # f = Future() - # f.set_result(self._calc_M(scratch, unsuccessful_solve_msg, upper)) - # running[constraint, other_disjunct, upper] = f + # f.set_result(self._calc_M(scratch, unsuccessful_solve_msg, is_upper)) + # running[constraint, other_disjunct, is_upper] = f to_deactivate = set() - for (constraint, other_disjunct, upper), future in running.items(): + for (constraint, other_disjunct, is_upper), future in running.items(): # try: M, disjunct_infeasible = future.result() # print(f"Got result: {M}, {disjunct_infeasible}") @@ -472,7 +487,7 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp continue # Note that we can't just transform immediately because we might be # waiting on the other one of upper_M or lower_M. - if upper: + if is_upper: Ms[constraint, other_disjunct] = (Ms[constraint, other_disjunct][0], M) else: Ms[constraint, other_disjunct] = (M, Ms[constraint, other_disjunct][1]) @@ -499,10 +514,10 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp disjunction._algebraic_constraint = weakref_ref(algebraic_constraint[disjunction.index()]) disjunction.deactivate() - def _get_scratch_block(self, constraint, other_disjunct, all_vars, upper): + def _get_scratch_block(self, constraint, other_disjunct, all_vars, is_upper): # Should this be a Block or a ConcreteModel? scratch = Block() - if upper: + if is_upper: scratch.obj = Objective(expr=constraint.body - constraint.upper, sense=maximize) else: scratch.obj = Objective(expr=constraint.body - constraint.lower, sense=minimize) @@ -519,7 +534,7 @@ def _get_scratch_block(self, constraint, other_disjunct, all_vars, upper): scratch.construct() # should this go in _calc_M? return scratch - def _calc_M(self, scratch, unsuccessful_message, upper): + def _calc_M(self, scratch, unsuccessful_message, is_upper): # Retrieve a solver object from thread-local storage; if we don't have one, make # a new solver reasonably similar to the one passed via config. if hasattr(_thread_local, 'solver'): @@ -528,9 +543,10 @@ def _calc_M(self, scratch, unsuccessful_message, upper): solver_arg = self._config.solver solver = _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) # print("About to call solver") - results = solver.solve(scratch, tee=False, load_solutions=False) + results = solver.solve(scratch, tee=True, load_solutions=False, keepfiles=False) if results.solver.termination_condition is TerminationCondition.infeasible: + # print("infeasible") # [2/18/24]: TODO: After the solver rewrite is complete, we will not # need this check since we can actually determine from the # termination condition whether or not the solver proved @@ -560,10 +576,50 @@ def _calc_M(self, scratch, unsuccessful_message, upper): # We should use the dual bound here. The primal bound is mathematically # incorrect in the presence of error, but local solvers need it. - if upper != self._config.use_primal_bound: - M = results.problem.upper_bound + if not self._config.use_primal_bound: + if is_upper: + M = results.problem.upper_bound + else: + M = results.problem.lower_bound + elif solver.name not in _solvers_require_loading_solutions: + if is_upper: + M = results.problem.lower_bound + else: + M = results.problem.upper_bound else: - M = results.problem.lower_bound + # In this unfortunate case, we checked earlier and ensured there is only + # one worker thread. It's obviously not safe to call load_from without a + # lock, but then the call to solver.solve would also need to acquire the + # same lock, otherwise the other threads would read a model halfway + # through a load_solutions call. But doing that would defeat the purpose + # of multithreading. Cloning `scratch` is not enough to fix this problem: + # the decision variables live on the original model, so we would have to + # clone the entire thing, which is likely expensive enough to negate any + # benefit gained by multithreading. + + # TODO: could I evaluate the objective without touching the vars by using + # replace_expressions and the var map in the results object? If so, will + # it be ipopt-only or will it work on any solver with this problem? + + # Non-multithreaded solution + # scratch.solutions.load_from(results) + # M = value(scratch.obj) + + # Potentially ipopt-specific workaround (do we always populate + # results._smap in the same way?) + M = value( + EXPR.replace_expressions( + scratch.obj, + substitution_map={ + id: results.solution.variable[name]['Value'] + for id, name in results._smap.byObject.items() + if isinstance(results._smap.bySymbol[name], Var) + }, + descend_into_named_expressions=True, + remove_named_expressions=True + ) + ) + if not math.isfinite(M): raise GDP_Error( "Subproblem solved to optimality, but no finite dual bound was " @@ -573,9 +629,10 @@ def _calc_M(self, scratch, unsuccessful_message, upper): "aware that interior feasible points for this subproblem do " "not give valid values for M." if not self._config.use_primal_bound else - "No finite objective after optimal solve." + "No finite objective after optimal solve. Possibly the solver " + "does not save the objective value in the results, in which case " + "it may need to be added to `_solvers_require_loading_solutions`" ) - # print(f"{M=}") return (M, False) def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms): From cdb0103a3bd9cd46aa9563d494c60e18ed7d56f4 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Thu, 12 Jun 2025 09:35:42 -0600 Subject: [PATCH 07/26] Handle ipopt in a somewhat better way --- pyomo/gdp/plugins/multiple_bigm.py | 126 +++++++++++++---------------- 1 file changed, 58 insertions(+), 68 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 948f6b9d081..0cbc71ed9ef 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -40,7 +40,7 @@ Var, ConcreteModel, ) -from pyomo.core.base import Reference, TransformationFactory +from pyomo.core.base import Reference, TransformationFactory, VarData import pyomo.core.expr as EXPR from pyomo.core.util import target_list from pyomo.common.errors import DeveloperError @@ -78,16 +78,6 @@ 'highs', } -# Solvers that do not report objective bounds or values in the result and therefore -# require loading the solution. For now, this prevents the use of multithreading; see the -# comment in _calc_M for more information. - -# TODO: this no longer prevents multithreading unless I go back to the singlethreaded -# method -_solvers_require_loading_solutions = { - 'ipopt', -} - # Keep Solver objects in thread-local storage _thread_local = threading.local() @@ -442,11 +432,6 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp Ms = arg_Ms # (Now exiting the DisjunctionDatas loop) threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) - ## Not doing this, unless I decide to do it again - # solver_name = self._config.solver.name - # if threads > 1 and solver_name in _solvers_require_loading_solutions: - # logger.warning(f"Solver '{solver_name}' cannot be used in multithreaded M calculations for now - reverting to a single thread. Pass the option threads=1 to silence this warning. ") - # threads = 1 running = {} if jobs: logger.info(f"Executing {len(jobs)} jobs on {threads} worker threads.") @@ -543,7 +528,7 @@ def _calc_M(self, scratch, unsuccessful_message, is_upper): solver_arg = self._config.solver solver = _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) # print("About to call solver") - results = solver.solve(scratch, tee=True, load_solutions=False, keepfiles=False) + results = solver.solve(scratch, tee=False, load_solutions=False, keepfiles=False) if results.solver.termination_condition is TerminationCondition.infeasible: # print("infeasible") @@ -571,68 +556,73 @@ def _calc_M(self, scratch, unsuccessful_message, is_upper): raise GDP_Error(unsuccessful_message) else: # print("We were optimal; getting M from objective bound") - # note: This tranformation can be made faster by allowing the solver a gap. As - # long as we have a bound, it's still valid (though not as tight). + + # note: This tranformation can be made faster by allowing the solver a + # gap. As long as we have a bound, it's still valid (though not as tight). # We should use the dual bound here. The primal bound is mathematically - # incorrect in the presence of error, but local solvers need it. + # incorrect in the presence of numerical error, but it's the best a local + # solver like ipopt can do, so we allow it to be used by setting an option. if not self._config.use_primal_bound: if is_upper: M = results.problem.upper_bound else: M = results.problem.lower_bound - elif solver.name not in _solvers_require_loading_solutions: + if not math.isfinite(M): + raise GDP_Error( + "Subproblem solved to optimality, but no finite dual bound was " + "obtained. If you would like to instead use the best obtained " + "solution, set the option use_primal_bound=True. This is " + "necessary when using a local solver such as ipopt, but be " + "aware that interior feasible points for this subproblem do " + "not give valid values for M in general." + ) + else: if is_upper: M = results.problem.lower_bound else: M = results.problem.upper_bound - else: - # In this unfortunate case, we checked earlier and ensured there is only - # one worker thread. It's obviously not safe to call load_from without a - # lock, but then the call to solver.solve would also need to acquire the - # same lock, otherwise the other threads would read a model halfway - # through a load_solutions call. But doing that would defeat the purpose - # of multithreading. Cloning `scratch` is not enough to fix this problem: - # the decision variables live on the original model, so we would have to - # clone the entire thing, which is likely expensive enough to negate any - # benefit gained by multithreading. - - # TODO: could I evaluate the objective without touching the vars by using - # replace_expressions and the var map in the results object? If so, will - # it be ipopt-only or will it work on any solver with this problem? - - # Non-multithreaded solution - # scratch.solutions.load_from(results) - # M = value(scratch.obj) - - # Potentially ipopt-specific workaround (do we always populate - # results._smap in the same way?) - M = value( - EXPR.replace_expressions( - scratch.obj, - substitution_map={ - id: results.solution.variable[name]['Value'] - for id, name in results._smap.byObject.items() - if isinstance(results._smap.bySymbol[name], Var) - }, - descend_into_named_expressions=True, - remove_named_expressions=True - ) - ) - - if not math.isfinite(M): - raise GDP_Error( - "Subproblem solved to optimality, but no finite dual bound was " - "obtained. If you would like to instead use the best obtained " - "solution, set the option use_primal_bound=True. This is " - "necessary when using a local solver such as ipopt, but be " - "aware that interior feasible points for this subproblem do " - "not give valid values for M." - if not self._config.use_primal_bound else - "No finite objective after optimal solve. Possibly the solver " - "does not save the objective value in the results, in which case " - "it may need to be added to `_solvers_require_loading_solutions`" - ) + if not math.isfinite(M): + # We are looking for the incumbent objective value after claiming to + # have solved the problem to optimality, but the primal bound on the + # results object is infinite or nan. It's possible something has gone + # very wrong, but it's more likely that we are using a solver + # interface, like ipopt, that rudely neglects to actually fill in + # that field for us. We will try again to get the value of the + # objective, assuming the results object has an _smap attribute that + # behaves like the one in pyomo.opt.base.solvers.OptSolver (as of + # 06/2025), and give up otherwise. + # + # The most robust way to do this would be to actually load the + # solution and only then evaluate the objective. But this isn't + # possible in the multithreaded context (solve() would read values of + # variables on the main model while we're writing), so we would have + # to switch to passing back the result object and waiting to parse + # until we're singlethreaded again, and I don't want to do that just + # to handle this edge case better. + try: + M = value( + EXPR.replace_expressions( + scratch.obj, + substitution_map={ + id: results.solution.variable[name]['Value'] + for id, name in results._smap.byObject.items() + if isinstance(results._smap.bySymbol[name], VarData) + }, + descend_into_named_expressions=True, + remove_named_expressions=True + ) + ) + if not math.isfinite(M): + raise ValueError() + except Exception: + raise GDP_Error( + "`use_primal_bound` is enabled, but could not find a finite" + "objective value after optimal solve. Possibly the solver " + "interface does not store the objective value in an easily " + "acessible location." + ) + # print(f"{M=}") return (M, False) def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms): From 160fc22493a6d6b9e928e6e2339152b2b7635a16 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Thu, 12 Jun 2025 09:48:52 -0600 Subject: [PATCH 08/26] apply black --- pyomo/gdp/plugins/multiple_bigm.py | 162 ++++++++++++++++++++--------- 1 file changed, 115 insertions(+), 47 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 0cbc71ed9ef..2931439b952 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -81,6 +81,7 @@ # Keep Solver objects in thread-local storage _thread_local = threading.local() + def Solver(val): if isinstance(val, str): return SolverFactory(val) @@ -235,8 +236,8 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): doc=""" If not specified, use up to the number of available cores, minus one (but at least one). - """ - ) + """, + ), ) CONFIG.declare( 'use_primal_bound', @@ -249,8 +250,8 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): aware that interior feasible points for this subproblem do not give valid values for M. That is, in the presence of any numerical error, this option will lead to a slightly wrong reformulation. - """ - ) + """, + ), ) transformation_name = 'mbigm' @@ -307,13 +308,17 @@ def _apply_to_impl(self, instance, **kwds): preprocessed_targets = gdp_tree.reverse_topological_sort() arg_Ms = self._config.bigM if self._config.bigM is not None else {} - self._transform_disjunctionDatas_threaded(preprocessed_targets, arg_Ms, gdp_tree) + self._transform_disjunctionDatas_threaded( + preprocessed_targets, arg_Ms, gdp_tree + ) # issue warnings about anything that was in the bigM args dict that we # didn't use _warn_for_unused_bigM_args(self._config.bigM, self.used_args, logger) - def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp_tree): + def _transform_disjunctionDatas_threaded( + self, preprocessed_targets, arg_Ms, gdp_tree + ): # We wish we could do this one Disjunction at a time, but we also want to # calculate the Ms in parallel. So we do a slightly convoluted dance of iterating # the Disjunctions to get a list of M calculation jobs, then calculating Ms @@ -361,35 +366,39 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp ) # start doing transformation - (transBlock, algebraic_constraint) = self._setup_transform_disjunctionData( - t, gdp_tree.root_disjunct(t) + (transBlock, algebraic_constraint) = ( + self._setup_transform_disjunctionData(t, gdp_tree.root_disjunct(t)) ) disjunction_setup[t] = (transBlock, algebraic_constraint) # already got arg_Ms - active_disjuncts[t] = ComponentSet(disj for disj in t.disjuncts if disj.active) + active_disjuncts[t] = ComponentSet( + disj for disj in t.disjuncts if disj.active + ) # this method returns the constraints transformed on this disjunct, # use update because we are now saving these from all disjuncts for # later if self._config.reduce_bound_constraints: - transformed_constraints.update(self._transform_bound_constraints( - active_disjuncts[t], transBlock, arg_Ms - )) + transformed_constraints.update( + self._transform_bound_constraints( + active_disjuncts[t], transBlock, arg_Ms + ) + ) # Get the jobs to calculate missing M values for this Disjunction. We can # skip this if we are only doing bound constraints. if not self._config.only_mbigm_bound_constraints: all_vars = list(self._get_all_var_objects(active_disjuncts[t])) for disjunct, other_disjunct in itertools.product( - active_disjuncts[t], active_disjuncts[t] + active_disjuncts[t], active_disjuncts[t] ): if disjunct is other_disjunct: continue for constraint in disjunct.component_data_objects( - Constraint, - active=True, - descend_into=Block, - sort=SortComponents.deterministic, + Constraint, + active=True, + descend_into=Block, + sort=SortComponents.deterministic, ): if constraint in transformed_constraints: continue @@ -397,9 +406,14 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp if (constraint, other_disjunct) in arg_Ms: # print("Getting Ms from arg") (lower_M, upper_M) = _convert_M_to_tuple( - arg_Ms[constraint, other_disjunct], constraint, other_disjunct + arg_Ms[constraint, other_disjunct], + constraint, + other_disjunct, + ) + self.used_args[constraint, other_disjunct] = ( + lower_M, + upper_M, ) - self.used_args[constraint, other_disjunct] = (lower_M, upper_M) else: # print("Did not find Ms in arg") (lower_M, upper_M) = (None, None) @@ -411,15 +425,38 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp ) if constraint.lower is not None and lower_M is None: - scratch = self._get_scratch_block(constraint, other_disjunct, all_vars, False) - jobs.append((constraint, other_disjunct, scratch, unsuccessful_solve_msg, False)) + scratch = self._get_scratch_block( + constraint, other_disjunct, all_vars, False + ) + jobs.append( + ( + constraint, + other_disjunct, + scratch, + unsuccessful_solve_msg, + False, + ) + ) if constraint.upper is not None and upper_M is None: - scratch = self._get_scratch_block(constraint, other_disjunct, all_vars, True) - jobs.append((constraint, other_disjunct, scratch, unsuccessful_solve_msg, True)) + scratch = self._get_scratch_block( + constraint, other_disjunct, all_vars, True + ) + jobs.append( + ( + constraint, + other_disjunct, + scratch, + unsuccessful_solve_msg, + True, + ) + ) # Let's avoid touching arg_Ms and just have one source of # truth for the current Ms Ms[constraint, other_disjunct] = (lower_M, upper_M) - transBlock._mbm_values[constraint, other_disjunct] = (lower_M, upper_M) + transBlock._mbm_values[constraint, other_disjunct] = ( + lower_M, + upper_M, + ) # (Now exiting: if not only_mbigm_bound_constraints) # # In this case we just pass through arg_Ms as-is and let @@ -431,19 +468,26 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp else: Ms = arg_Ms # (Now exiting the DisjunctionDatas loop) - threads = self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) + threads = ( + self._config.threads + if self._config.threads is not None + else max(1, len(os.sched_getaffinity(0)) - 1) + ) running = {} if jobs: logger.info(f"Executing {len(jobs)} jobs on {threads} worker threads.") - with ThreadPoolExecutor(max_workers=threads, thread_name_prefix="mBigM-Worker-Thread") as executor: - for constraint, other_disjunct, scratch, unsuccessful_solve_msg, is_upper in jobs: - running[constraint, other_disjunct, is_upper] = ( - executor.submit( - self._calc_M, - scratch, - unsuccessful_solve_msg, - is_upper, - ) + with ThreadPoolExecutor( + max_workers=threads, thread_name_prefix="mBigM-Worker-Thread" + ) as executor: + for ( + constraint, + other_disjunct, + scratch, + unsuccessful_solve_msg, + is_upper, + ) in jobs: + running[constraint, other_disjunct, is_upper] = executor.submit( + self._calc_M, scratch, unsuccessful_solve_msg, is_upper ) # # cheat for testing purposes # f = Future() @@ -466,17 +510,26 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp # If we made it here without an exception, the solver is on the # trusted solvers list logger.debug( - "Disjunct '%s' is infeasible, deactivating." % other_disjunct.name + "Disjunct '%s' is infeasible, deactivating." + % other_disjunct.name ) to_deactivate.add(other_disjunct) continue # Note that we can't just transform immediately because we might be # waiting on the other one of upper_M or lower_M. if is_upper: - Ms[constraint, other_disjunct] = (Ms[constraint, other_disjunct][0], M) + Ms[constraint, other_disjunct] = ( + Ms[constraint, other_disjunct][0], + M, + ) else: - Ms[constraint, other_disjunct] = (M, Ms[constraint, other_disjunct][1]) - transBlock._mbm_values[constraint, other_disjunct] = Ms[constraint, other_disjunct] + Ms[constraint, other_disjunct] = ( + M, + Ms[constraint, other_disjunct][1], + ) + transBlock._mbm_values[constraint, other_disjunct] = Ms[ + constraint, other_disjunct + ] # Now exiting threaded context. Cleanup: for disjunct in to_deactivate: @@ -488,24 +541,35 @@ def _transform_disjunctionDatas_threaded(self, preprocessed_targets, arg_Ms, gdp # print(f"Final {Ms.values()=}") # Iterate the Disjunctions again to actually transform them - for disjunction, (transBlock, algebraic_constraint) in disjunction_setup.items(): + for disjunction, ( + transBlock, + algebraic_constraint, + ) in disjunction_setup.items(): or_expr = 0 for disjunct in active_disjuncts[disjunction]: or_expr += disjunct.indicator_var.get_associated_binary() - self._transform_disjunct(disjunct, transBlock, active_disjuncts[disjunction], Ms) + self._transform_disjunct( + disjunct, transBlock, active_disjuncts[disjunction], Ms + ) algebraic_constraint.add(disjunction.index(), or_expr == 1) # map the DisjunctionData to its XOR constraint to mark it as # transformed - disjunction._algebraic_constraint = weakref_ref(algebraic_constraint[disjunction.index()]) + disjunction._algebraic_constraint = weakref_ref( + algebraic_constraint[disjunction.index()] + ) disjunction.deactivate() - + def _get_scratch_block(self, constraint, other_disjunct, all_vars, is_upper): # Should this be a Block or a ConcreteModel? scratch = Block() if is_upper: - scratch.obj = Objective(expr=constraint.body - constraint.upper, sense=maximize) + scratch.obj = Objective( + expr=constraint.body - constraint.upper, sense=maximize + ) else: - scratch.obj = Objective(expr=constraint.body - constraint.lower, sense=minimize) + scratch.obj = Objective( + expr=constraint.body - constraint.lower, sense=minimize + ) # Create a Block component (reference) that actually points to a DisjunctData, as # we want the writer to write it as if it were an ordinary Block rather than # getting any special handling. DisjunctData inherits BlockData, so this should @@ -526,9 +590,13 @@ def _calc_M(self, scratch, unsuccessful_message, is_upper): solver = _thread_local.solver else: solver_arg = self._config.solver - solver = _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) + solver = _thread_local.solver = SolverFactory( + solver_arg.name, options=solver_arg.options + ) # print("About to call solver") - results = solver.solve(scratch, tee=False, load_solutions=False, keepfiles=False) + results = solver.solve( + scratch, tee=False, load_solutions=False, keepfiles=False + ) if results.solver.termination_condition is TerminationCondition.infeasible: # print("infeasible") @@ -610,7 +678,7 @@ def _calc_M(self, scratch, unsuccessful_message, is_upper): if isinstance(results._smap.bySymbol[name], VarData) }, descend_into_named_expressions=True, - remove_named_expressions=True + remove_named_expressions=True, ) ) if not math.isfinite(M): From 83a2a20fa42f8623ea4d7f6871c31e8d3dc2ad79 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Thu, 12 Jun 2025 10:04:31 -0600 Subject: [PATCH 09/26] apply black on more files --- pyomo/gdp/plugins/bigm_mixin.py | 1 + pyomo/gdp/tests/test_mbigm.py | 5 ++++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/bigm_mixin.py b/pyomo/gdp/plugins/bigm_mixin.py index a8cd9051ed4..dc9a1cbf42d 100644 --- a/pyomo/gdp/plugins/bigm_mixin.py +++ b/pyomo/gdp/plugins/bigm_mixin.py @@ -19,6 +19,7 @@ logger = logging.getLogger('pyomo.gdp.bigm') + def _convert_M_to_tuple(M, constraint, disjunct=None): if not isinstance(M, (tuple, list)): if M is None: diff --git a/pyomo/gdp/tests/test_mbigm.py b/pyomo/gdp/tests/test_mbigm.py index 56daf21e77f..34b6f9e4eb2 100644 --- a/pyomo/gdp/tests/test_mbigm.py +++ b/pyomo/gdp/tests/test_mbigm.py @@ -1062,7 +1062,10 @@ def test_calculate_Ms_infeasible_Disjunct_local_solver(self): r"Disjunct 'disjunction_disjuncts\[0\]' is selected.", ): TransformationFactory('gdp.mbigm').apply_to( - m, solver=SolverFactory('ipopt'), reduce_bound_constraints=False, use_primal_bound=True + m, + solver=SolverFactory('ipopt'), + reduce_bound_constraints=False, + use_primal_bound=True, ) @unittest.skipUnless(gurobi_available, "Gurobi is not available") From 2e5e07e975d17b113049f29952cecbdd5e578c92 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Thu, 12 Jun 2025 12:25:55 -0600 Subject: [PATCH 10/26] Cleanup and formatting --- pyomo/gdp/plugins/multiple_bigm.py | 318 ++++++++++++++--------------- 1 file changed, 156 insertions(+), 162 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 2931439b952..88eea4c90c0 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -61,7 +61,7 @@ from weakref import ref as weakref_ref import threading -from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor, Future +from concurrent.futures import ThreadPoolExecutor import os logger = logging.getLogger('pyomo.gdp.mbigm') @@ -308,17 +308,13 @@ def _apply_to_impl(self, instance, **kwds): preprocessed_targets = gdp_tree.reverse_topological_sort() arg_Ms = self._config.bigM if self._config.bigM is not None else {} - self._transform_disjunctionDatas_threaded( - preprocessed_targets, arg_Ms, gdp_tree - ) + self._transform_disjunctionDatas(preprocessed_targets, arg_Ms, gdp_tree) # issue warnings about anything that was in the bigM args dict that we # didn't use _warn_for_unused_bigM_args(self._config.bigM, self.used_args, logger) - def _transform_disjunctionDatas_threaded( - self, preprocessed_targets, arg_Ms, gdp_tree - ): + def _transform_disjunctionDatas(self, preprocessed_targets, arg_Ms, gdp_tree): # We wish we could do this one Disjunction at a time, but we also want to # calculate the Ms in parallel. So we do a slightly convoluted dance of iterating # the Disjunctions to get a list of M calculation jobs, then calculating Ms @@ -366,18 +362,16 @@ def _transform_disjunctionDatas_threaded( ) # start doing transformation - (transBlock, algebraic_constraint) = ( + disjunction_setup[t] = (transBlock, algebraic_constraint) = ( self._setup_transform_disjunctionData(t, gdp_tree.root_disjunct(t)) ) - disjunction_setup[t] = (transBlock, algebraic_constraint) # already got arg_Ms active_disjuncts[t] = ComponentSet( disj for disj in t.disjuncts if disj.active ) - # this method returns the constraints transformed on this disjunct, - # use update because we are now saving these from all disjuncts for - # later + # this method returns the constraints transformed on this disjunct; + # update because we are saving these from all disjuncts for later if self._config.reduce_bound_constraints: transformed_constraints.update( self._transform_bound_constraints( @@ -385,151 +379,77 @@ def _transform_disjunctionDatas_threaded( ) ) # Get the jobs to calculate missing M values for this Disjunction. We can - # skip this if we are only doing bound constraints. + # skip this if we are only doing bound constraints, in which case we pass + # through arg_Ms as-is and let transform_constraint parse it. if not self._config.only_mbigm_bound_constraints: - all_vars = list(self._get_all_var_objects(active_disjuncts[t])) - for disjunct, other_disjunct in itertools.product( - active_disjuncts[t], active_disjuncts[t] - ): - if disjunct is other_disjunct: - continue - - for constraint in disjunct.component_data_objects( - Constraint, - active=True, - descend_into=Block, - sort=SortComponents.deterministic, - ): - if constraint in transformed_constraints: - continue - # First check args - if (constraint, other_disjunct) in arg_Ms: - # print("Getting Ms from arg") - (lower_M, upper_M) = _convert_M_to_tuple( - arg_Ms[constraint, other_disjunct], - constraint, - other_disjunct, - ) - self.used_args[constraint, other_disjunct] = ( - lower_M, - upper_M, - ) - else: - # print("Did not find Ms in arg") - (lower_M, upper_M) = (None, None) - unsuccessful_solve_msg = ( - "Unsuccessful solve to calculate M value to " - "relax constraint '%s' on Disjunct '%s' when " - "Disjunct '%s' is selected." - % (constraint.name, disjunct.name, other_disjunct.name) - ) - - if constraint.lower is not None and lower_M is None: - scratch = self._get_scratch_block( - constraint, other_disjunct, all_vars, False - ) - jobs.append( - ( - constraint, - other_disjunct, - scratch, - unsuccessful_solve_msg, - False, - ) - ) - if constraint.upper is not None and upper_M is None: - scratch = self._get_scratch_block( - constraint, other_disjunct, all_vars, True - ) - jobs.append( - ( - constraint, - other_disjunct, - scratch, - unsuccessful_solve_msg, - True, - ) - ) - # Let's avoid touching arg_Ms and just have one source of - # truth for the current Ms - Ms[constraint, other_disjunct] = (lower_M, upper_M) - transBlock._mbm_values[constraint, other_disjunct] = ( - lower_M, - upper_M, - ) - # (Now exiting: if not only_mbigm_bound_constraints) - # - # In this case we just pass through arg_Ms as-is and let - # transform_constraint parse it - # - # TODO revisit what's going on here. Can we just skip half the - # transformation instead of doing this? We have queued zero jobs so - # none of the threads stuff will run. + self._setup_jobs_for_disjunction( + t, + active_disjuncts, + transformed_constraints, + arg_Ms, + Ms, + jobs, + transBlock, + ) else: Ms = arg_Ms # (Now exiting the DisjunctionDatas loop) - threads = ( - self._config.threads - if self._config.threads is not None - else max(1, len(os.sched_getaffinity(0)) - 1) - ) - running = {} + to_deactivate = set() if jobs: - logger.info(f"Executing {len(jobs)} jobs on {threads} worker threads.") - with ThreadPoolExecutor( - max_workers=threads, thread_name_prefix="mBigM-Worker-Thread" - ) as executor: - for ( - constraint, - other_disjunct, - scratch, - unsuccessful_solve_msg, - is_upper, - ) in jobs: - running[constraint, other_disjunct, is_upper] = executor.submit( - self._calc_M, scratch, unsuccessful_solve_msg, is_upper - ) - # # cheat for testing purposes - # f = Future() - # f.set_result(self._calc_M(scratch, unsuccessful_solve_msg, is_upper)) - # running[constraint, other_disjunct, is_upper] = f - - to_deactivate = set() - for (constraint, other_disjunct, is_upper), future in running.items(): - # try: - M, disjunct_infeasible = future.result() - # print(f"Got result: {M}, {disjunct_infeasible}") - # except Exception as e: - # import traceback - # print("@@@@@@@@@@@ PRINTING ERROR! @@@@@@@@@") - # traceback.print_tb(e.__traceback__) - # raise e - if disjunct_infeasible: - # Don't deactivate live because we have concurrent jobs still running - if other_disjunct not in to_deactivate: - # If we made it here without an exception, the solver is on the - # trusted solvers list - logger.debug( - "Disjunct '%s' is infeasible, deactivating." - % other_disjunct.name - ) - to_deactivate.add(other_disjunct) - continue - # Note that we can't just transform immediately because we might be - # waiting on the other one of upper_M or lower_M. - if is_upper: - Ms[constraint, other_disjunct] = ( - Ms[constraint, other_disjunct][0], - M, - ) - else: - Ms[constraint, other_disjunct] = ( - M, - Ms[constraint, other_disjunct][1], + threads = ( + self._config.threads + if self._config.threads is not None + else max(1, len(os.sched_getaffinity(0)) - 1) + ) + running = {} + logger.info(f"Running {len(jobs)} jobs on {threads} worker threads.") + with ThreadPoolExecutor( + max_workers=threads, thread_name_prefix="mBigM-Worker-Thread" + ) as executor: + for ( + constraint, + other_disjunct, + scratch, + unsuccessful_solve_msg, + is_upper, + ) in jobs: + running[constraint, other_disjunct, is_upper] = executor.submit( + self._calc_M, scratch, unsuccessful_solve_msg, is_upper ) - transBlock._mbm_values[constraint, other_disjunct] = Ms[ - constraint, other_disjunct - ] + # # cheat for testing purposes + # f = Future() + # f.set_result(self._calc_M(scratch, unsuccessful_solve_msg, is_upper)) + # running[constraint, other_disjunct, is_upper] = f + + for (constraint, other_disjunct, is_upper), future in running.items(): + M, disjunct_infeasible = future.result() + if disjunct_infeasible: + # Don't deactivate live because we have concurrent jobs still + # running + if other_disjunct not in to_deactivate: + # If we made it here without an exception, the solver is on + # the trusted solvers list + logger.debug( + "Disjunct '%s' is infeasible, deactivating." + % other_disjunct.name + ) + to_deactivate.add(other_disjunct) + continue + # Note that we can't just transform immediately because we might be + # waiting on the other one of upper_M or lower_M. + if is_upper: + Ms[constraint, other_disjunct] = ( + Ms[constraint, other_disjunct][0], + M, + ) + else: + Ms[constraint, other_disjunct] = ( + M, + Ms[constraint, other_disjunct][1], + ) + transBlock._mbm_values[constraint, other_disjunct] = Ms[ + constraint, other_disjunct + ] # Now exiting threaded context. Cleanup: for disjunct in to_deactivate: @@ -539,7 +459,6 @@ def _transform_disjunctionDatas_threaded( for con in transformed_constraints: con.deactivate() - # print(f"Final {Ms.values()=}") # Iterate the Disjunctions again to actually transform them for disjunction, ( transBlock, @@ -583,23 +502,102 @@ def _get_scratch_block(self, constraint, other_disjunct, all_vars, is_upper): scratch.construct() # should this go in _calc_M? return scratch + # Do the inner work setting up M calculation jobs for one disjunction. This mutates + # the parameters Ms, jobs, transBlock._mbm_values, and also self.used_args. + def _setup_jobs_for_disjunction( + self, + disjunction, + active_disjuncts, + transformed_constraints, + arg_Ms, + Ms, + jobs, + transBlock, + ): + all_vars = list(self._get_all_var_objects(active_disjuncts[disjunction])) + for disjunct, other_disjunct in itertools.product( + active_disjuncts[disjunction], active_disjuncts[disjunction] + ): + if disjunct is other_disjunct: + continue + + for constraint in disjunct.component_data_objects( + Constraint, + active=True, + descend_into=Block, + sort=SortComponents.deterministic, + ): + if constraint in transformed_constraints: + continue + # First check args + if (constraint, other_disjunct) in arg_Ms: + (lower_M, upper_M) = _convert_M_to_tuple( + arg_Ms[constraint, other_disjunct], + constraint, + other_disjunct, + ) + self.used_args[constraint, other_disjunct] = ( + lower_M, + upper_M, + ) + else: + (lower_M, upper_M) = (None, None) + unsuccessful_solve_msg = ( + "Unsuccessful solve to calculate M value to " + "relax constraint '%s' on Disjunct '%s' when " + "Disjunct '%s' is selected." + % (constraint.name, disjunct.name, other_disjunct.name) + ) + + if constraint.lower is not None and lower_M is None: + scratch = self._get_scratch_block( + constraint, other_disjunct, all_vars, False + ) + jobs.append( + ( + constraint, + other_disjunct, + scratch, + unsuccessful_solve_msg, + False, + ) + ) + if constraint.upper is not None and upper_M is None: + scratch = self._get_scratch_block( + constraint, other_disjunct, all_vars, True + ) + jobs.append( + ( + constraint, + other_disjunct, + scratch, + unsuccessful_solve_msg, + True, + ) + ) + + Ms[constraint, other_disjunct] = (lower_M, upper_M) + transBlock._mbm_values[constraint, other_disjunct] = ( + lower_M, + upper_M, + ) + def _calc_M(self, scratch, unsuccessful_message, is_upper): # Retrieve a solver object from thread-local storage; if we don't have one, make - # a new solver reasonably similar to the one passed via config. - if hasattr(_thread_local, 'solver'): + # a new solver reasonably similar to the one passed via config. We never actually + # use the one we were passed. + if hasattr(_thread_local, "solver"): solver = _thread_local.solver else: solver_arg = self._config.solver solver = _thread_local.solver = SolverFactory( solver_arg.name, options=solver_arg.options ) - # print("About to call solver") results = solver.solve( scratch, tee=False, load_solutions=False, keepfiles=False ) if results.solver.termination_condition is TerminationCondition.infeasible: - # print("infeasible") # [2/18/24]: TODO: After the solver rewrite is complete, we will not # need this check since we can actually determine from the # termination condition whether or not the solver proved @@ -620,12 +618,9 @@ def _calc_M(self, scratch, unsuccessful_message, is_upper): # and then we can abandon this hack raise GDP_Error(unsuccessful_message) elif results.solver.termination_condition is not TerminationCondition.optimal: - # print("We had a solver error!") raise GDP_Error(unsuccessful_message) else: - # print("We were optimal; getting M from objective bound") - - # note: This tranformation can be made faster by allowing the solver a + # note: This transformation can be made faster by allowing the solver a # gap. As long as we have a bound, it's still valid (though not as tight). # We should use the dual bound here. The primal bound is mathematically @@ -688,9 +683,8 @@ def _calc_M(self, scratch, unsuccessful_message, is_upper): "`use_primal_bound` is enabled, but could not find a finite" "objective value after optimal solve. Possibly the solver " "interface does not store the objective value in an easily " - "acessible location." + "accessible location." ) - # print(f"{M=}") return (M, False) def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms): From 21ea319647ad596f905ed9c0abae467be4ca4561 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Thu, 12 Jun 2025 14:16:59 -0600 Subject: [PATCH 11/26] black yet again --- pyomo/gdp/plugins/multiple_bigm.py | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 88eea4c90c0..0df3b406b84 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -532,14 +532,9 @@ def _setup_jobs_for_disjunction( # First check args if (constraint, other_disjunct) in arg_Ms: (lower_M, upper_M) = _convert_M_to_tuple( - arg_Ms[constraint, other_disjunct], - constraint, - other_disjunct, - ) - self.used_args[constraint, other_disjunct] = ( - lower_M, - upper_M, + arg_Ms[constraint, other_disjunct], constraint, other_disjunct ) + self.used_args[constraint, other_disjunct] = (lower_M, upper_M) else: (lower_M, upper_M) = (None, None) unsuccessful_solve_msg = ( @@ -577,10 +572,7 @@ def _setup_jobs_for_disjunction( ) Ms[constraint, other_disjunct] = (lower_M, upper_M) - transBlock._mbm_values[constraint, other_disjunct] = ( - lower_M, - upper_M, - ) + transBlock._mbm_values[constraint, other_disjunct] = (lower_M, upper_M) def _calc_M(self, scratch, unsuccessful_message, is_upper): # Retrieve a solver object from thread-local storage; if we don't have one, make From c9794a5ab70708d12d1d72a53005831725657aff Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 18 Jun 2025 11:29:23 -0600 Subject: [PATCH 12/26] Switch from multithreading to multiprocessing --- pyomo/gdp/plugins/multiple_bigm.py | 514 ++++++++++++++++------------- 1 file changed, 292 insertions(+), 222 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 0df3b406b84..8d9146551b1 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -40,7 +40,7 @@ Var, ConcreteModel, ) -from pyomo.core.base import Reference, TransformationFactory, VarData +from pyomo.core.base import Reference, TransformationFactory import pyomo.core.expr as EXPR from pyomo.core.util import target_list from pyomo.common.errors import DeveloperError @@ -60,9 +60,10 @@ from weakref import ref as weakref_ref -import threading -from concurrent.futures import ThreadPoolExecutor +import multiprocessing import os +import threading +from pyomo.common.dependencies import dill, dill_available logger = logging.getLogger('pyomo.gdp.mbigm') @@ -78,8 +79,14 @@ 'highs', } -# Keep Solver objects in thread-local storage +# Whether these are thread-local or at module scope makes no difference for 'spawn' or +# 'forkserver', but it matters if we run single-threaded, and possibly for edge cases of +# 'fork' (although it's not correct to use fork() here while having multiple threads +# anyway, making it moot in theory). _thread_local = threading.local() +_thread_local.solver = None +_thread_local.model = None +_thread_local.config_use_primal_bound = None def Solver(val): @@ -90,6 +97,13 @@ def Solver(val): return val +def StartMethod(val): + if val in {'spawn', 'fork', 'forkserver', None}: + return val + else: + raise ValueError("Expected one of 'spawn', 'fork', or 'forkserver', or None.") + + @TransformationFactory.register( 'gdp.mbigm', doc="Relax disjunctive model using big-M terms specific to each disjunct", @@ -145,7 +159,7 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): CONFIG.declare( 'solver', ConfigValue( - default='baron', # TODO restore this to gurobi after fixing deadlock + default='gurobi', domain=Solver, description="A solver to use to solve the continuous subproblems for " "calculating the M values", @@ -232,10 +246,29 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): ConfigValue( default=None, domain=int, - description="Number of worker threads to use when estimating M values", + description="Number of worker processes to use when estimating M values", doc=""" If not specified, use up to the number of available cores, minus - one (but at least one). + one (but at least one). If set to zero, revert to single-threaded operation. + """, + ), + ) + CONFIG.declare( + 'process_start_method', + ConfigValue( + default=None, + domain=StartMethod, + description="Start method used for spawning processes during M calculation", + doc=""" + Options are 'fork', 'spawn', and 'forkserver'. See the Python + multiprocessing documentation for a full description of each of these. When + None is passed, we determine a reasonable default. On POSIX, the default is + 'fork', unless we detect that Python has multiple threads at the time the + process pool is created, in which case we instead use 'forkserver'. On + Windows, the default and only possible option is 'spawn'. Note that if + 'spawn' or 'forkserver' are selected, we depend on the `dill` module for + pickling, and model instances must be pickleable using `dill`. This option + is ignored if `threads` is set to 0. """, ), ) @@ -244,11 +277,12 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): ConfigValue( default=False, domain=bool, - description="When estimating M values, use the primal bound instead of dual bound", + description="When estimating M values, use the primal bound " + "instead of the dual bound.", doc=""" This is necessary when using a local solver such as ipopt, but be aware that interior feasible points for this subproblem do not give - valid values for M. That is, in the presence of any numerical error, + valid values for M. That is, in the presence of numerical error, this option will lead to a slightly wrong reformulation. """, ), @@ -272,6 +306,9 @@ def _apply_to(self, instance, **kwds): self._arg_list.clear() self._expr_bound_visitor.leaf_bounds.clear() self._expr_bound_visitor.use_fixed_var_values_as_bounds = False + _thread_local.model = None + _thread_local.solver = None + _thread_local.config_use_primal_bound = None def _apply_to_impl(self, instance, **kwds): self._process_arguments(instance, **kwds) @@ -308,28 +345,35 @@ def _apply_to_impl(self, instance, **kwds): preprocessed_targets = gdp_tree.reverse_topological_sort() arg_Ms = self._config.bigM if self._config.bigM is not None else {} - self._transform_disjunctionDatas(preprocessed_targets, arg_Ms, gdp_tree) + self._transform_disjunctionDatas( + instance, preprocessed_targets, arg_Ms, gdp_tree + ) # issue warnings about anything that was in the bigM args dict that we # didn't use _warn_for_unused_bigM_args(self._config.bigM, self.used_args, logger) - def _transform_disjunctionDatas(self, preprocessed_targets, arg_Ms, gdp_tree): + def _transform_disjunctionDatas( + self, instance, preprocessed_targets, arg_Ms, gdp_tree + ): # We wish we could do this one Disjunction at a time, but we also want to # calculate the Ms in parallel. So we do a slightly convoluted dance of iterating # the Disjunctions to get a list of M calculation jobs, then calculating Ms # multithreaded, then iterating Disjunctions again to actually transform the # constraints. - # To-do list in form (constraint, other_disj, scratch, unsucc_msg, is_upper) + # To-do list in form (constraint, other_disjunct, unsuccessful_message, is_upper) jobs = [] # map Disjunction -> set of its active Disjuncts active_disjuncts = ComponentMap() # set of Constraints processed during special handling of bound constraints: we # will deactivate these, but not until we're done calculating Ms transformed_constraints = set() - # Finished M values + # Finished M values. If we are only doing the bound constraints, we will skip all + # the calculation steps and pass these through to transform_constraints. Ms = {} + if self._config.only_mbigm_bound_constraints: + Ms = arg_Ms # Disjunctions and their setup components. We will return to these after # calculating Ms. disjunction_setup = {} @@ -366,7 +410,6 @@ def _transform_disjunctionDatas(self, preprocessed_targets, arg_Ms, gdp_tree): self._setup_transform_disjunctionData(t, gdp_tree.root_disjunct(t)) ) - # already got arg_Ms active_disjuncts[t] = ComponentSet( disj for disj in t.disjuncts if disj.active ) @@ -378,9 +421,9 @@ def _transform_disjunctionDatas(self, preprocessed_targets, arg_Ms, gdp_tree): active_disjuncts[t], transBlock, arg_Ms ) ) - # Get the jobs to calculate missing M values for this Disjunction. We can - # skip this if we are only doing bound constraints, in which case we pass - # through arg_Ms as-is and let transform_constraint parse it. + # Get the jobs to calculate missing M values for this Disjunction. We + # skip this if we are only doing bound constraints, in which case Ms was + # already set. if not self._config.only_mbigm_bound_constraints: self._setup_jobs_for_disjunction( t, @@ -391,75 +434,125 @@ def _transform_disjunctionDatas(self, preprocessed_targets, arg_Ms, gdp_tree): jobs, transBlock, ) - else: - Ms = arg_Ms # (Now exiting the DisjunctionDatas loop) - to_deactivate = set() if jobs: threads = ( self._config.threads if self._config.threads is not None else max(1, len(os.sched_getaffinity(0)) - 1) ) - running = {} - logger.info(f"Running {len(jobs)} jobs on {threads} worker threads.") - with ThreadPoolExecutor( - max_workers=threads, thread_name_prefix="mBigM-Worker-Thread" - ) as executor: - for ( - constraint, - other_disjunct, - scratch, - unsuccessful_solve_msg, - is_upper, - ) in jobs: - running[constraint, other_disjunct, is_upper] = executor.submit( - self._calc_M, scratch, unsuccessful_solve_msg, is_upper + if threads > 0: + method = ( + self._config.process_start_method + if self._config.process_start_method is not None + else ( + 'spawn' + if os.name == 'nt' + else 'forkserver' if len(threading.enumerate()) > 1 else 'fork' ) - # # cheat for testing purposes - # f = Future() - # f.set_result(self._calc_M(scratch, unsuccessful_solve_msg, is_upper)) - # running[constraint, other_disjunct, is_upper] = f - - for (constraint, other_disjunct, is_upper), future in running.items(): - M, disjunct_infeasible = future.result() - if disjunct_infeasible: - # Don't deactivate live because we have concurrent jobs still - # running - if other_disjunct not in to_deactivate: - # If we made it here without an exception, the solver is on - # the trusted solvers list - logger.debug( - "Disjunct '%s' is infeasible, deactivating." - % other_disjunct.name - ) - to_deactivate.add(other_disjunct) - continue - # Note that we can't just transform immediately because we might be - # waiting on the other one of upper_M or lower_M. - if is_upper: - Ms[constraint, other_disjunct] = ( - Ms[constraint, other_disjunct][0], - M, + ) + logger.info( + f"Running {len(jobs)} jobs on {threads} " + f"worker processes with method {method}." + ) + if method == 'spawn' or method == 'forkserver': + if not dill_available: + raise GDP_Error( + "Dill is required when spawning processes using " + "methods 'spawn' or 'forkserver', but it could " + "not be imported." ) - else: - Ms[constraint, other_disjunct] = ( - M, - Ms[constraint, other_disjunct][1], + pool = multiprocessing.get_context(method).Pool( + processes=threads, + initializer=_setup_spawn, + initargs=( + dill.dumps(instance), + dill.dumps(self._config.solver), + self._config.use_primal_bound, + ), + ) + elif method == 'fork': + _thread_local.model = instance + _thread_local.solver = self._config.solver + _thread_local.config_use_primal_bound = ( + self._config.use_primal_bound + ) + pool = multiprocessing.get_context('fork').Pool( + processes=threads, initializer=_setup_fork, initargs=() + ) + with pool: + results = pool.starmap( + func=_calc_M, + iterable=[ + ( + constraint.getname(fully_qualified=True), + other_disjunct.getname(fully_qualified=True), + unsuccessful_solve_msg, + is_upper, + ) + for ( + constraint, + other_disjunct, + unsuccessful_solve_msg, + is_upper, + ) in jobs + ], + ) + else: + _thread_local.model = instance + _thread_local.solver = self._config.solver + results = itertools.starmap( + _calc_M, + [ + ( + constraint.getname(fully_qualified=True), + other_disjunct.getname(fully_qualified=True), + unsuccessful_solve_msg, + is_upper, ) - transBlock._mbm_values[constraint, other_disjunct] = Ms[ - constraint, other_disjunct - ] - - # Now exiting threaded context. Cleanup: - for disjunct in to_deactivate: - disjunct.deactivate() - active_disjuncts[gdp_tree.parent(disjunct)].remove(disjunct) + for ( + constraint, + other_disjunct, + unsuccessful_solve_msg, + is_upper, + ) in jobs + ], + ) + for (constraint, other_disjunct, _, is_upper), ( + M, + disjunct_infeasible, + ) in zip(jobs, results): + if disjunct_infeasible: + # If we made it here without an exception, the solver is on the + # trusted solvers list + logger.debug( + "Disjunct '%s' is infeasible, deactivating." + % other_disjunct.name + ) + other_disjunct.deactivate() + active_disjuncts[gdp_tree.parent(other_disjunct)].remove( + other_disjunct + ) + # Note that we can't just transform immediately because we might be + # waiting on the other one of upper_M or lower_M. + if is_upper: + Ms[constraint, other_disjunct] = ( + Ms[constraint, other_disjunct][0], + M, + ) + else: + Ms[constraint, other_disjunct] = ( + M, + Ms[constraint, other_disjunct][1], + ) + transBlock._mbm_values[constraint, other_disjunct] = Ms[ + constraint, other_disjunct + ] for con in transformed_constraints: con.deactivate() - # Iterate the Disjunctions again to actually transform them + # Iterate the Disjunctions again and actually transform them for disjunction, ( transBlock, algebraic_constraint, @@ -478,30 +571,6 @@ def _transform_disjunctionDatas(self, preprocessed_targets, arg_Ms, gdp_tree): ) disjunction.deactivate() - def _get_scratch_block(self, constraint, other_disjunct, all_vars, is_upper): - # Should this be a Block or a ConcreteModel? - scratch = Block() - if is_upper: - scratch.obj = Objective( - expr=constraint.body - constraint.upper, sense=maximize - ) - else: - scratch.obj = Objective( - expr=constraint.body - constraint.lower, sense=minimize - ) - # Create a Block component (reference) that actually points to a DisjunctData, as - # we want the writer to write it as if it were an ordinary Block rather than - # getting any special handling. DisjunctData inherits BlockData, so this should - # be valid. - scratch.disjunct_ref = Reference(other_disjunct, ctype=Block) - # If the writers don't assume Vars are declared on the Block - # being solved, we won't need this! - for v in all_vars: - ref = Reference(v) - scratch.add_component(unique_component_name(scratch, v.name), ref) - scratch.construct() # should this go in _calc_M? - return scratch - # Do the inner work setting up M calculation jobs for one disjunction. This mutates # the parameters Ms, jobs, transBlock._mbm_values, and also self.used_args. def _setup_jobs_for_disjunction( @@ -514,7 +583,6 @@ def _setup_jobs_for_disjunction( jobs, transBlock, ): - all_vars = list(self._get_all_var_objects(active_disjuncts[disjunction])) for disjunct, other_disjunct in itertools.product( active_disjuncts[disjunction], active_disjuncts[disjunction] ): @@ -545,140 +613,17 @@ def _setup_jobs_for_disjunction( ) if constraint.lower is not None and lower_M is None: - scratch = self._get_scratch_block( - constraint, other_disjunct, all_vars, False - ) jobs.append( - ( - constraint, - other_disjunct, - scratch, - unsuccessful_solve_msg, - False, - ) + (constraint, other_disjunct, unsuccessful_solve_msg, False) ) if constraint.upper is not None and upper_M is None: - scratch = self._get_scratch_block( - constraint, other_disjunct, all_vars, True - ) jobs.append( - ( - constraint, - other_disjunct, - scratch, - unsuccessful_solve_msg, - True, - ) + (constraint, other_disjunct, unsuccessful_solve_msg, True) ) Ms[constraint, other_disjunct] = (lower_M, upper_M) transBlock._mbm_values[constraint, other_disjunct] = (lower_M, upper_M) - def _calc_M(self, scratch, unsuccessful_message, is_upper): - # Retrieve a solver object from thread-local storage; if we don't have one, make - # a new solver reasonably similar to the one passed via config. We never actually - # use the one we were passed. - if hasattr(_thread_local, "solver"): - solver = _thread_local.solver - else: - solver_arg = self._config.solver - solver = _thread_local.solver = SolverFactory( - solver_arg.name, options=solver_arg.options - ) - results = solver.solve( - scratch, tee=False, load_solutions=False, keepfiles=False - ) - - if results.solver.termination_condition is TerminationCondition.infeasible: - # [2/18/24]: TODO: After the solver rewrite is complete, we will not - # need this check since we can actually determine from the - # termination condition whether or not the solver proved - # infeasibility or just terminated at local infeasiblity. For now, - # while this is not complete, it catches most of the solvers we - # trust, and, unless someone is so pathological as to *rename* an - # untrusted solver using a trusted solver name, it will never do the - # *wrong* thing. - if any(s in solver.name for s in _trusted_solvers): - return (0, True) - else: - # This is a solver that might report - # 'infeasible' for local infeasibility, so we - # can't deactivate with confidence. To be - # conservative, we'll just complain about - # it. Post-solver-rewrite we will want to change - # this so that we check for 'proven_infeasible' - # and then we can abandon this hack - raise GDP_Error(unsuccessful_message) - elif results.solver.termination_condition is not TerminationCondition.optimal: - raise GDP_Error(unsuccessful_message) - else: - # note: This transformation can be made faster by allowing the solver a - # gap. As long as we have a bound, it's still valid (though not as tight). - - # We should use the dual bound here. The primal bound is mathematically - # incorrect in the presence of numerical error, but it's the best a local - # solver like ipopt can do, so we allow it to be used by setting an option. - if not self._config.use_primal_bound: - if is_upper: - M = results.problem.upper_bound - else: - M = results.problem.lower_bound - if not math.isfinite(M): - raise GDP_Error( - "Subproblem solved to optimality, but no finite dual bound was " - "obtained. If you would like to instead use the best obtained " - "solution, set the option use_primal_bound=True. This is " - "necessary when using a local solver such as ipopt, but be " - "aware that interior feasible points for this subproblem do " - "not give valid values for M in general." - ) - else: - if is_upper: - M = results.problem.lower_bound - else: - M = results.problem.upper_bound - if not math.isfinite(M): - # We are looking for the incumbent objective value after claiming to - # have solved the problem to optimality, but the primal bound on the - # results object is infinite or nan. It's possible something has gone - # very wrong, but it's more likely that we are using a solver - # interface, like ipopt, that rudely neglects to actually fill in - # that field for us. We will try again to get the value of the - # objective, assuming the results object has an _smap attribute that - # behaves like the one in pyomo.opt.base.solvers.OptSolver (as of - # 06/2025), and give up otherwise. - # - # The most robust way to do this would be to actually load the - # solution and only then evaluate the objective. But this isn't - # possible in the multithreaded context (solve() would read values of - # variables on the main model while we're writing), so we would have - # to switch to passing back the result object and waiting to parse - # until we're singlethreaded again, and I don't want to do that just - # to handle this edge case better. - try: - M = value( - EXPR.replace_expressions( - scratch.obj, - substitution_map={ - id: results.solution.variable[name]['Value'] - for id, name in results._smap.byObject.items() - if isinstance(results._smap.bySymbol[name], VarData) - }, - descend_into_named_expressions=True, - remove_named_expressions=True, - ) - ) - if not math.isfinite(M): - raise ValueError() - except Exception: - raise GDP_Error( - "`use_primal_bound` is enabled, but could not find a finite" - "objective value after optimal solve. Possibly the solver " - "interface does not store the objective value in an easily " - "accessible location." - ) - return (M, False) - def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms): # We've already filtered out deactivated disjuncts, so we know obj is # active. @@ -972,3 +917,128 @@ def get_all_M_values(self, model): all_ms.update(transBlock._mbm_values) return all_ms + + +# Things we call in subprocesses. These can't be member functions, or else we'd have to +# pickle `self`, which is problematic. +def _setup_spawn(model, solver, use_primal_bound): + _thread_local.model = dill.loads(model) + solver_arg = dill.loads(solver) + # I'm not sure if this is necessary after it's been pickled, but let's be on the safe + # side + # _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) + _thread_local.solver = solver_arg + _thread_local.config_use_primal_bound = use_primal_bound + + +def _setup_fork(): + # model and config_use_primal_bound were already properly set, but remake the solver + # instead of using the passed argument. All these processes are copies of the calling + # thread so the thread-local still works. + _thread_local.solver = SolverFactory( + _thread_local.solver.name, options=_thread_local.solver.options + ) + + +def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper): + solver = _thread_local.solver + model = _thread_local.model + scratch = _get_scratch_block( + model.find_component(constraint_name), + model.find_component(other_disjunct_name), + is_upper, + ) + results = solver.solve(scratch, tee=False, load_solutions=False, keepfiles=False) + + if results.solver.termination_condition is TerminationCondition.infeasible: + # [2/18/24]: TODO: After the solver rewrite is complete, we will not + # need this check since we can actually determine from the + # termination condition whether or not the solver proved + # infeasibility or just terminated at local infeasiblity. For now, + # while this is not complete, it catches most of the solvers we + # trust, and, unless someone is so pathological as to *rename* an + # untrusted solver using a trusted solver name, it will never do the + # *wrong* thing. + if any(s in solver.name for s in _trusted_solvers): + return (0, True) + else: + # This is a solver that might report + # 'infeasible' for local infeasibility, so we + # can't deactivate with confidence. To be + # conservative, we'll just complain about + # it. Post-solver-rewrite we will want to change + # this so that we check for 'proven_infeasible' + # and then we can abandon this hack + raise GDP_Error(unsuccessful_message) + elif results.solver.termination_condition is not TerminationCondition.optimal: + raise GDP_Error(unsuccessful_message) + else: + # note: This transformation can be made faster by allowing the solver a + # gap. As long as we have a bound, it's still valid (though not as tight). + + # We should use the dual bound here. The primal bound is mathematically + # incorrect in the presence of numerical error, but it's the best a local + # solver like ipopt can do, so we allow it to be used by setting an option. + if not _thread_local.config_use_primal_bound: + if is_upper: + M = results.problem.upper_bound + else: + M = results.problem.lower_bound + if not math.isfinite(M): + raise GDP_Error( + "Subproblem solved to optimality, but no finite dual bound was " + "obtained. If you would like to instead use the best obtained " + "solution, set the option use_primal_bound=True. This is " + "necessary when using a local solver such as ipopt, but be " + "aware that interior feasible points for this subproblem do " + "not give valid values for M in general." + ) + else: + if is_upper: + M = results.problem.lower_bound + else: + M = results.problem.upper_bound + if not math.isfinite(M): + # Solved to optimality, but we did find an incumbent objective value. Try + # again by actually loading the solution and evaluating the objective + # expression. + try: + scratch.solutions.load_from(results) + M = value(scratch.obj) + if not math.isfinite(M): + raise ValueError() + except Exception: + raise GDP_Error( + "`use_primal_bound` is enabled, but could not find a finite" + "objective value after optimal solve." + ) + return (M, False) + + +def _get_scratch_block(constraint, other_disjunct, is_upper): + # Should this be a Block or a ConcreteModel? + scratch = Block() + if is_upper: + scratch.obj = Objective(expr=constraint.body - constraint.upper, sense=maximize) + else: + scratch.obj = Objective(expr=constraint.body - constraint.lower, sense=minimize) + # Create a Block component (reference) that actually points to a DisjunctData, as + # we want the writer to write it as if it were an ordinary Block rather than + # getting any special handling. DisjunctData inherits BlockData, so this should + # be valid. + scratch.disjunct_ref = Reference(other_disjunct, ctype=Block) + + # Add references to every Var that appears in an active constraint. + # TODO: If the writers don't assume Vars are declared on the Block + # being solved, we won't need this! + seen = set() + for constraint in scratch.component_data_objects( + Constraint, active=True, sort=SortComponents.deterministic, descend_into=Block + ): + for var in EXPR.identify_variables(constraint.expr, include_fixed=True): + if id(var) not in seen: + seen.add(id(var)) + ref = Reference(var) + scratch.add_component(unique_component_name(scratch, var.name), ref) + scratch.construct() + return scratch From a2d52cffc2873d80220da67099eb51671b593a96 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 18 Jun 2025 12:12:48 -0600 Subject: [PATCH 13/26] fix a bug, also reflow my block comments --- pyomo/gdp/plugins/multiple_bigm.py | 103 ++++++++++++++++------------- 1 file changed, 57 insertions(+), 46 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 8d9146551b1..7a8d896642a 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -79,10 +79,11 @@ 'highs', } -# Whether these are thread-local or at module scope makes no difference for 'spawn' or -# 'forkserver', but it matters if we run single-threaded, and possibly for edge cases of -# 'fork' (although it's not correct to use fork() here while having multiple threads -# anyway, making it moot in theory). +# Whether these are thread-local or at module scope makes no difference +# for 'spawn' or 'forkserver', but it matters if we run single-threaded, +# and possibly for edge cases of 'fork' (although it's not correct to +# use fork() here while having multiple threads anyway, making it moot +# in theory). _thread_local = threading.local() _thread_local.solver = None _thread_local.model = None @@ -356,21 +357,24 @@ def _apply_to_impl(self, instance, **kwds): def _transform_disjunctionDatas( self, instance, preprocessed_targets, arg_Ms, gdp_tree ): - # We wish we could do this one Disjunction at a time, but we also want to - # calculate the Ms in parallel. So we do a slightly convoluted dance of iterating - # the Disjunctions to get a list of M calculation jobs, then calculating Ms - # multithreaded, then iterating Disjunctions again to actually transform the + # We wish we could do this one Disjunction at a time, but we + # also want to calculate the Ms in parallel. So we do a slightly + # convoluted dance of iterating the Disjunctions to get a list + # of M calculation jobs, then calculating Ms multithreaded, then + # iterating Disjunctions again to actually transform the # constraints. # To-do list in form (constraint, other_disjunct, unsuccessful_message, is_upper) jobs = [] # map Disjunction -> set of its active Disjuncts active_disjuncts = ComponentMap() - # set of Constraints processed during special handling of bound constraints: we - # will deactivate these, but not until we're done calculating Ms + # set of Constraints processed during special handling of bound + # constraints: we will deactivate these, but not until we're + # done calculating Ms transformed_constraints = set() - # Finished M values. If we are only doing the bound constraints, we will skip all - # the calculation steps and pass these through to transform_constraints. + # Finished M values. If we are only doing the bound constraints, + # we will skip all the calculation steps and pass these through + # to transform_constraints. Ms = {} if self._config.only_mbigm_bound_constraints: Ms = arg_Ms @@ -518,21 +522,24 @@ def _transform_disjunctionDatas( ) in jobs ], ) + deactivated = set() for (constraint, other_disjunct, _, is_upper), ( M, disjunct_infeasible, ) in zip(jobs, results): if disjunct_infeasible: - # If we made it here without an exception, the solver is on the - # trusted solvers list - logger.debug( - "Disjunct '%s' is infeasible, deactivating." - % other_disjunct.name - ) - other_disjunct.deactivate() - active_disjuncts[gdp_tree.parent(other_disjunct)].remove( - other_disjunct - ) + if other_disjunct not in deactivated: + # If we made it here without an exception, the solver is on the + # trusted solvers list + logger.debug( + "Disjunct '%s' is infeasible, deactivating." + % other_disjunct.name + ) + other_disjunct.deactivate() + active_disjuncts[gdp_tree.parent(other_disjunct)].remove( + other_disjunct + ) + deactivated.add(other_disjunct) # Note that we can't just transform immediately because we might be # waiting on the other one of upper_M or lower_M. if is_upper: @@ -571,8 +578,9 @@ def _transform_disjunctionDatas( ) disjunction.deactivate() - # Do the inner work setting up M calculation jobs for one disjunction. This mutates - # the parameters Ms, jobs, transBlock._mbm_values, and also self.used_args. + # Do the inner work setting up M calculation jobs for one + # disjunction. This mutates the parameters Ms, jobs, + # transBlock._mbm_values, and also self.used_args. def _setup_jobs_for_disjunction( self, disjunction, @@ -919,22 +927,22 @@ def get_all_M_values(self, model): return all_ms -# Things we call in subprocesses. These can't be member functions, or else we'd have to -# pickle `self`, which is problematic. +# Things we call in subprocesses. These can't be member functions, or +# else we'd have to pickle `self`, which is problematic. def _setup_spawn(model, solver, use_primal_bound): _thread_local.model = dill.loads(model) solver_arg = dill.loads(solver) - # I'm not sure if this is necessary after it's been pickled, but let's be on the safe - # side - # _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) - _thread_local.solver = solver_arg + # I'm not sure if this is necessary after it's been pickled, but + # let's be on the safe side + _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) _thread_local.config_use_primal_bound = use_primal_bound def _setup_fork(): - # model and config_use_primal_bound were already properly set, but remake the solver - # instead of using the passed argument. All these processes are copies of the calling - # thread so the thread-local still works. + # model and config_use_primal_bound were already properly set, but + # remake the solver instead of using the passed argument. All these + # processes are copies of the calling thread so the thread-local + # still works. _thread_local.solver = SolverFactory( _thread_local.solver.name, options=_thread_local.solver.options ) @@ -973,12 +981,14 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper elif results.solver.termination_condition is not TerminationCondition.optimal: raise GDP_Error(unsuccessful_message) else: - # note: This transformation can be made faster by allowing the solver a - # gap. As long as we have a bound, it's still valid (though not as tight). - - # We should use the dual bound here. The primal bound is mathematically - # incorrect in the presence of numerical error, but it's the best a local - # solver like ipopt can do, so we allow it to be used by setting an option. + # NOTE: This transformation can be made faster by allowing the + # solver a gap. As long as we have a bound, it's still valid + # (though not as tight). + # + # We should use the dual bound here. The primal bound is + # mathematically incorrect in the presence of numerical error, + # but it's the best a local solver like ipopt can do, so we + # allow it to be used by setting an option. if not _thread_local.config_use_primal_bound: if is_upper: M = results.problem.upper_bound @@ -999,9 +1009,9 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper else: M = results.problem.upper_bound if not math.isfinite(M): - # Solved to optimality, but we did find an incumbent objective value. Try - # again by actually loading the solution and evaluating the objective - # expression. + # Solved to optimality, but we did not find an incumbent + # objective value. Try again by actually loading the + # solution and evaluating the objective expression. try: scratch.solutions.load_from(results) M = value(scratch.obj) @@ -1022,10 +1032,11 @@ def _get_scratch_block(constraint, other_disjunct, is_upper): scratch.obj = Objective(expr=constraint.body - constraint.upper, sense=maximize) else: scratch.obj = Objective(expr=constraint.body - constraint.lower, sense=minimize) - # Create a Block component (reference) that actually points to a DisjunctData, as - # we want the writer to write it as if it were an ordinary Block rather than - # getting any special handling. DisjunctData inherits BlockData, so this should - # be valid. + # Create a Block component (reference) that actually points to a + # DisjunctData, as we want the writer to write it as if it were an + # ordinary Block rather than getting any special + # handling. DisjunctData inherits BlockData, so this should be + # valid. scratch.disjunct_ref = Reference(other_disjunct, ctype=Block) # Add references to every Var that appears in an active constraint. From 738223b40a23422c091aa5b254bc9a4da7279611 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 18 Jun 2025 14:44:36 -0600 Subject: [PATCH 14/26] Import plugins manually, also remove unused imports --- pyomo/gdp/plugins/multiple_bigm.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 7a8d896642a..2fd6d422396 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -19,31 +19,20 @@ from pyomo.common.modeling import unique_component_name from pyomo.core import ( - Any, - Binary, Block, - BooleanVar, - Connector, Constraint, - Expression, - ExternalFunction, maximize, minimize, NonNegativeIntegers, Objective, - Param, - Set, - SetOf, SortComponents, Suffix, value, - Var, - ConcreteModel, + Any, ) from pyomo.core.base import Reference, TransformationFactory import pyomo.core.expr as EXPR from pyomo.core.util import target_list -from pyomo.common.errors import DeveloperError from pyomo.gdp import Disjunct, Disjunction, GDP_Error from pyomo.gdp.plugins.bigm_mixin import ( @@ -52,9 +41,7 @@ _warn_for_unused_bigM_args, ) from pyomo.gdp.plugins.gdp_to_mip_transformation import GDP_to_MIP_Transformation -from pyomo.gdp.transformed_disjunct import _TransformedDisjunct -from pyomo.gdp.util import get_gdp_tree, _to_dict -from pyomo.network import Port +from pyomo.gdp.util import _to_dict from pyomo.opt import SolverFactory, TerminationCondition from pyomo.repn import generate_standard_repn @@ -65,6 +52,18 @@ import threading from pyomo.common.dependencies import dill, dill_available +# When using 'spawn' or 'forkserver', Python starts in a new environment +# and executes only this file, so we need to manually ensure necessary +# plugins are registered (even if the calling process has already +# registered them). +import pyomo.solvers.plugins +import pyomo.opt.plugins +import pyomo.repn.plugins + +pyomo.solvers.plugins.load() +pyomo.opt.plugins.load() +pyomo.repn.plugins.load() + logger = logging.getLogger('pyomo.gdp.mbigm') _trusted_solvers = { @@ -455,6 +454,7 @@ def _transform_disjunctionDatas( else 'forkserver' if len(threading.enumerate()) > 1 else 'fork' ) ) + method = 'spawn' logger.info( f"Running {len(jobs)} jobs on {threads} " f"worker processes with method {method}." From bb7a816a988196c546148f5d748d6ef037af053b Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 18 Jun 2025 14:46:01 -0600 Subject: [PATCH 15/26] remove debug code --- pyomo/gdp/plugins/multiple_bigm.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 2fd6d422396..ad218a54f52 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -454,7 +454,6 @@ def _transform_disjunctionDatas( else 'forkserver' if len(threading.enumerate()) > 1 else 'fork' ) ) - method = 'spawn' logger.info( f"Running {len(jobs)} jobs on {threads} " f"worker processes with method {method}." From 1c0bb3ecadcd8336a2add331db7658186a6cda48 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 20 Jun 2025 13:35:34 -0600 Subject: [PATCH 16/26] Set 'name' class attribute of LegacySolver objects to the legacy name --- pyomo/contrib/solver/common/factory.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyomo/contrib/solver/common/factory.py b/pyomo/contrib/solver/common/factory.py index 29aca32b7cb..0320062a925 100644 --- a/pyomo/contrib/solver/common/factory.py +++ b/pyomo/contrib/solver/common/factory.py @@ -31,6 +31,7 @@ def decorator(cls): class LegacySolver(LegacySolverWrapper, cls): pass + LegacySolver.name = legacy_name LegacySolverFactory.register(legacy_name, doc + " (new interface)")( LegacySolver ) From d57501d3620faaf82e96d24e5863369e5165b298 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 20 Jun 2025 14:44:03 -0600 Subject: [PATCH 17/26] Consider and fix more edge cases --- pyomo/gdp/plugins/multiple_bigm.py | 97 +++++++++++++++--------------- pyomo/gdp/tests/test_mbigm.py | 9 ++- 2 files changed, 56 insertions(+), 50 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index ad218a54f52..e65b7b82ad1 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -20,6 +20,7 @@ from pyomo.core import ( Block, + ConcreteModel, Constraint, maximize, minimize, @@ -43,6 +44,8 @@ from pyomo.gdp.plugins.gdp_to_mip_transformation import GDP_to_MIP_Transformation from pyomo.gdp.util import _to_dict from pyomo.opt import SolverFactory, TerminationCondition +from pyomo.contrib.solver.common.base import SolverBase as NewSolverBase +from pyomo.contrib.solver.common.base import LegacySolverWrapper from pyomo.repn import generate_standard_repn from weakref import ref as weakref_ref @@ -52,17 +55,6 @@ import threading from pyomo.common.dependencies import dill, dill_available -# When using 'spawn' or 'forkserver', Python starts in a new environment -# and executes only this file, so we need to manually ensure necessary -# plugins are registered (even if the calling process has already -# registered them). -import pyomo.solvers.plugins -import pyomo.opt.plugins -import pyomo.repn.plugins - -pyomo.solvers.plugins.load() -pyomo.opt.plugins.load() -pyomo.repn.plugins.load() logger = logging.getLogger('pyomo.gdp.mbigm') @@ -94,10 +86,15 @@ def Solver(val): return SolverFactory(val) if not hasattr(val, 'solve'): raise ValueError("Expected a string or solver object (with solve() method)") + if isinstance(val, NewSolverBase) and not isinstance(val, LegacySolverWrapper): + raise ValueError( + "Please pass an old-style solver object, using the " + "LegacySolverWrapper mechanism if necessary." + ) return val -def StartMethod(val): +def ProcessStartMethod(val): if val in {'spawn', 'fork', 'forkserver', None}: return val else: @@ -249,7 +246,8 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): description="Number of worker processes to use when estimating M values", doc=""" If not specified, use up to the number of available cores, minus - one (but at least one). If set to zero, revert to single-threaded operation. + one. If set to one or less, do not spawn processes, and revert to + single-threaded operation. """, ), ) @@ -257,10 +255,10 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): 'process_start_method', ConfigValue( default=None, - domain=StartMethod, + domain=ProcessStartMethod, description="Start method used for spawning processes during M calculation", doc=""" - Options are 'fork', 'spawn', and 'forkserver'. See the Python + Options are 'fork', 'spawn', and 'forkserver', or None. See the Python multiprocessing documentation for a full description of each of these. When None is passed, we determine a reasonable default. On POSIX, the default is 'fork', unless we detect that Python has multiple threads at the time the @@ -268,7 +266,7 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): Windows, the default and only possible option is 'spawn'. Note that if 'spawn' or 'forkserver' are selected, we depend on the `dill` module for pickling, and model instances must be pickleable using `dill`. This option - is ignored if `threads` is set to 0. + is ignored if `threads` is set to one or less. """, ), ) @@ -357,11 +355,11 @@ def _transform_disjunctionDatas( self, instance, preprocessed_targets, arg_Ms, gdp_tree ): # We wish we could do this one Disjunction at a time, but we - # also want to calculate the Ms in parallel. So we do a slightly - # convoluted dance of iterating the Disjunctions to get a list - # of M calculation jobs, then calculating Ms multithreaded, then - # iterating Disjunctions again to actually transform the - # constraints. + # also want to calculate the Ms in parallel. So we first iterate + # the Disjunctions once to get a list of M calculation jobs, + # then we calculate the Ms in parallel, then we return to a + # single thread and iterate Disjunctions again to actually + # transform the constraints. # To-do list in form (constraint, other_disjunct, unsuccessful_message, is_upper) jobs = [] @@ -442,9 +440,11 @@ def _transform_disjunctionDatas( threads = ( self._config.threads if self._config.threads is not None - else max(1, len(os.sched_getaffinity(0)) - 1) + # It would be better to use len(os.sched_getaffinity(0)), + # but it is not available on all platforms. + else os.cpu_count() - 1 ) - if threads > 0: + if threads > 1: method = ( self._config.process_start_method if self._config.process_start_method is not None @@ -455,8 +455,8 @@ def _transform_disjunctionDatas( ) ) logger.info( - f"Running {len(jobs)} jobs on {threads} " - f"worker processes with method {method}." + f"Running {len(jobs)} jobs on {threads} worker " + f"processes with process start method {method}." ) if method == 'spawn' or method == 'forkserver': if not dill_available: @@ -470,7 +470,8 @@ def _transform_disjunctionDatas( initializer=_setup_spawn, initargs=( dill.dumps(instance), - dill.dumps(self._config.solver), + self._config.solver.name, + self._config.solver.options, self._config.use_primal_bound, ), ) @@ -504,6 +505,7 @@ def _transform_disjunctionDatas( else: _thread_local.model = instance _thread_local.solver = self._config.solver + _thread_local.config_use_primal_bound = self._config.use_primal_bound results = itertools.starmap( _calc_M, [ @@ -928,12 +930,15 @@ def get_all_M_values(self, model): # Things we call in subprocesses. These can't be member functions, or # else we'd have to pickle `self`, which is problematic. -def _setup_spawn(model, solver, use_primal_bound): +def _setup_spawn(model, solver_name, solver_options, use_primal_bound): + # When using 'spawn' or 'forkserver', Python starts in a new + # environment and executes only this file, so we need to manually + # ensure necessary plugins are registered (even if the main process + # has already registered them). + import pyomo.environ + _thread_local.model = dill.loads(model) - solver_arg = dill.loads(solver) - # I'm not sure if this is necessary after it's been pickled, but - # let's be on the safe side - _thread_local.solver = SolverFactory(solver_arg.name, options=solver_arg.options) + _thread_local.solver = SolverFactory(solver_name, options=solver_options) _thread_local.config_use_primal_bound = use_primal_bound @@ -956,8 +961,14 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper is_upper, ) results = solver.solve(scratch, tee=False, load_solutions=False, keepfiles=False) - - if results.solver.termination_condition is TerminationCondition.infeasible: + termination_condition = results.solver.termination_condition + if is_upper: + incumbent = results.problem.lower_bound + bound = results.problem.upper_bound + else: + incumbent = results.problem.upper_bound + bound = results.problem.lower_bound + if termination_condition is TerminationCondition.infeasible: # [2/18/24]: TODO: After the solver rewrite is complete, we will not # need this check since we can actually determine from the # termination condition whether or not the solver proved @@ -977,7 +988,7 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper # this so that we check for 'proven_infeasible' # and then we can abandon this hack raise GDP_Error(unsuccessful_message) - elif results.solver.termination_condition is not TerminationCondition.optimal: + elif termination_condition is not TerminationCondition.optimal: raise GDP_Error(unsuccessful_message) else: # NOTE: This transformation can be made faster by allowing the @@ -989,10 +1000,7 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper # but it's the best a local solver like ipopt can do, so we # allow it to be used by setting an option. if not _thread_local.config_use_primal_bound: - if is_upper: - M = results.problem.upper_bound - else: - M = results.problem.lower_bound + M = bound if not math.isfinite(M): raise GDP_Error( "Subproblem solved to optimality, but no finite dual bound was " @@ -1003,10 +1011,7 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper "not give valid values for M in general." ) else: - if is_upper: - M = results.problem.lower_bound - else: - M = results.problem.upper_bound + M = incumbent if not math.isfinite(M): # Solved to optimality, but we did not find an incumbent # objective value. Try again by actually loading the @@ -1018,20 +1023,19 @@ def _calc_M(constraint_name, other_disjunct_name, unsuccessful_message, is_upper raise ValueError() except Exception: raise GDP_Error( - "`use_primal_bound` is enabled, but could not find a finite" + "`use_primal_bound` is enabled, but could not find a finite " "objective value after optimal solve." ) return (M, False) def _get_scratch_block(constraint, other_disjunct, is_upper): - # Should this be a Block or a ConcreteModel? - scratch = Block() + scratch = ConcreteModel() if is_upper: scratch.obj = Objective(expr=constraint.body - constraint.upper, sense=maximize) else: scratch.obj = Objective(expr=constraint.body - constraint.lower, sense=minimize) - # Create a Block component (reference) that actually points to a + # Create a Block component (via Reference) that actually points to a # DisjunctData, as we want the writer to write it as if it were an # ordinary Block rather than getting any special # handling. DisjunctData inherits BlockData, so this should be @@ -1050,5 +1054,4 @@ def _get_scratch_block(constraint, other_disjunct, is_upper): seen.add(id(var)) ref = Reference(var) scratch.add_component(unique_component_name(scratch, var.name), ref) - scratch.construct() return scratch diff --git a/pyomo/gdp/tests/test_mbigm.py b/pyomo/gdp/tests/test_mbigm.py index 34b6f9e4eb2..71f461bb83a 100644 --- a/pyomo/gdp/tests/test_mbigm.py +++ b/pyomo/gdp/tests/test_mbigm.py @@ -1054,12 +1054,15 @@ def test_calculate_Ms_infeasible_Disjunct(self): ) def test_calculate_Ms_infeasible_Disjunct_local_solver(self): m = self.make_infeasible_disjunct_model() + # When multiple exceptions are raised during a + # multiprocessing.Pool.map call, it is indeterminate which + # exception will be raised to the caller. with self.assertRaisesRegex( GDP_Error, r"Unsuccessful solve to calculate M value to " - r"relax constraint 'disjunction_disjuncts\[1\].constraint\[1\]' " - r"on Disjunct 'disjunction_disjuncts\[1\]' when " - r"Disjunct 'disjunction_disjuncts\[0\]' is selected.", + r"relax constraint 'disjunction_disjuncts\[\d+\].constraint\[\d+\]' " + r"on Disjunct 'disjunction_disjuncts\[\d+\]' when " + r"Disjunct 'disjunction_disjuncts\[\d+\]' is selected.", ): TransformationFactory('gdp.mbigm').apply_to( m, From 2de93fbc3ac4e36989b18207a4ed4844a0f06dad Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 20 Jun 2025 15:30:04 -0600 Subject: [PATCH 18/26] Test process spawning methods --- pyomo/gdp/plugins/multiple_bigm.py | 17 +------ pyomo/gdp/tests/test_mbigm.py | 80 +++++++++++++++++++++++++++++- 2 files changed, 80 insertions(+), 17 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index e65b7b82ad1..84b2bd7201d 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -506,6 +506,7 @@ def _transform_disjunctionDatas( _thread_local.model = instance _thread_local.solver = self._config.solver _thread_local.config_use_primal_bound = self._config.use_primal_bound + logger.info(f"Running {len(jobs)} jobs single-threaded.") results = itertools.starmap( _calc_M, [ @@ -849,22 +850,6 @@ def _add_transformation_block(self, to_block): ) return transBlock, new_block - def _get_all_var_objects(self, active_disjuncts): - # This is actually a general utility for getting all Vars that appear in - # active Disjuncts in a Disjunction. - seen = set() - for disj in active_disjuncts: - for constraint in disj.component_data_objects( - Constraint, - active=True, - sort=SortComponents.deterministic, - descend_into=Block, - ): - for var in EXPR.identify_variables(constraint.expr, include_fixed=True): - if id(var) not in seen: - seen.add(id(var)) - yield var - def _warn_for_active_suffix(self, suffix, disjunct, active_disjuncts, Ms): if suffix.local_name == 'BigM': logger.debug( diff --git a/pyomo/gdp/tests/test_mbigm.py b/pyomo/gdp/tests/test_mbigm.py index 71f461bb83a..37e11251978 100644 --- a/pyomo/gdp/tests/test_mbigm.py +++ b/pyomo/gdp/tests/test_mbigm.py @@ -13,6 +13,7 @@ import logging from os.path import join, normpath import pickle +import os from pyomo.common.fileutils import import_file, PYOMO_ROOT_DIR from pyomo.common.log import LoggingIntercept @@ -44,6 +45,7 @@ from pyomo.gdp.tests.models import make_indexed_equality_model from pyomo.repn import generate_standard_repn + gurobi_available = ( SolverFactory('gurobi').available(exception_flag=False) and SolverFactory('gurobi').license_is_valid() @@ -334,7 +336,7 @@ def check_linear_func_constraints(self, m, mbm, Ms=None): @unittest.skipUnless(gurobi_available, "Gurobi is not available") def test_calculated_Ms_correct(self): # Calculating all the Ms is expensive, so we just do it in this one test - # and then specify them for the others + # and then specify them for most of the others m = self.make_model() mbm = TransformationFactory('gdp.mbigm') mbm.apply_to(m, reduce_bound_constraints=False) @@ -906,6 +908,82 @@ def test_only_multiple_bigm_bound_constraints_arg_Ms(self): m, mbm, {m.d1: (-1050, 1050), m.d2: (-2000, 1200), m.d3: (-4000, 4000)} ) + # A set of tests identical to test_calculated_Ms_correct, except + # that we use each possible process spawning method for + # multiprocessing + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_calculated_Ms_spawn(self): + m = self.make_model() + mbm = TransformationFactory('gdp.mbigm') + mbm.apply_to( + m, reduce_bound_constraints=False, threads=3, process_start_method='spawn' + ) + + self.check_all_untightened_bounds_constraints(m, mbm) + self.check_linear_func_constraints(m, mbm) + + self.assertStructuredAlmostEqual(mbm.get_all_M_values(m), self.get_Ms(m)) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_calculated_Ms_singlethreaded(self): + m = self.make_model() + mbm = TransformationFactory('gdp.mbigm') + mbm.apply_to(m, reduce_bound_constraints=False, threads=1) + + self.check_all_untightened_bounds_constraints(m, mbm) + self.check_linear_func_constraints(m, mbm) + + self.assertStructuredAlmostEqual(mbm.get_all_M_values(m), self.get_Ms(m)) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + @unittest.skipIf(os.name == 'nt', "'forkserver' is not available on Windows") + def test_calculated_Ms_forkserver(self): + m = self.make_model() + mbm = TransformationFactory('gdp.mbigm') + mbm.apply_to( + m, + reduce_bound_constraints=False, + threads=3, + process_start_method='forkserver', + ) + + self.check_all_untightened_bounds_constraints(m, mbm) + self.check_linear_func_constraints(m, mbm) + + self.assertStructuredAlmostEqual(mbm.get_all_M_values(m), self.get_Ms(m)) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + @unittest.skipIf(os.name == 'nt', "'fork' is not available on Windows") + def test_calculated_Ms_fork(self): + m = self.make_model() + mbm = TransformationFactory('gdp.mbigm') + mbm.apply_to( + m, reduce_bound_constraints=False, threads=3, process_start_method='fork' + ) + + self.check_all_untightened_bounds_constraints(m, mbm) + self.check_linear_func_constraints(m, mbm) + + self.assertStructuredAlmostEqual(mbm.get_all_M_values(m), self.get_Ms(m)) + + # Make sure we don't choke on a LegacySolverWrapper + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_calculated_Ms_legacy_solver_wrapper(self): + m = self.make_model() + mbm = TransformationFactory('gdp.mbigm') + mbm.apply_to( + m, + reduce_bound_constraints=False, + threads=3, + process_start_method='spawn', + solver=SolverFactory('gurobi_direct_v2'), + ) + + self.check_all_untightened_bounds_constraints(m, mbm) + self.check_linear_func_constraints(m, mbm) + + self.assertStructuredAlmostEqual(mbm.get_all_M_values(m), self.get_Ms(m)) + @unittest.skipUnless(gurobi_available, "Gurobi is not available") class NestedDisjunctsInFlatGDP(unittest.TestCase): From 2c5d3e5be23fe273b25921e2da8b3f8a1d2e921f Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 20 Jun 2025 15:43:49 -0600 Subject: [PATCH 19/26] reorder some imports --- pyomo/gdp/plugins/multiple_bigm.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 84b2bd7201d..49334a07870 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -12,11 +12,15 @@ import itertools import logging import math +import multiprocessing +import os +import threading from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.common.config import ConfigDict, ConfigValue from pyomo.common.gc_manager import PauseGC from pyomo.common.modeling import unique_component_name +from pyomo.common.dependencies import dill, dill_available from pyomo.core import ( Block, @@ -50,11 +54,6 @@ from weakref import ref as weakref_ref -import multiprocessing -import os -import threading -from pyomo.common.dependencies import dill, dill_available - logger = logging.getLogger('pyomo.gdp.mbigm') From bcbb9ab8587ccbe5095cc2ceafcce03a49c8c354 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Jun 2025 13:59:24 -0600 Subject: [PATCH 20/26] Revert "Set 'name' class attribute of LegacySolver objects to the legacy name" This reverts commit 1c0bb3ecadcd8336a2add331db7658186a6cda48. --- pyomo/contrib/solver/common/factory.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyomo/contrib/solver/common/factory.py b/pyomo/contrib/solver/common/factory.py index 0320062a925..29aca32b7cb 100644 --- a/pyomo/contrib/solver/common/factory.py +++ b/pyomo/contrib/solver/common/factory.py @@ -31,7 +31,6 @@ def decorator(cls): class LegacySolver(LegacySolverWrapper, cls): pass - LegacySolver.name = legacy_name LegacySolverFactory.register(legacy_name, doc + " (new interface)")( LegacySolver ) From edc0f306e8d7d8bbe7f848baeb7057dad74e2e4c Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Jun 2025 14:19:55 -0600 Subject: [PATCH 21/26] Pickle classes instead of passing solver names --- pyomo/gdp/plugins/multiple_bigm.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 49334a07870..f1a8138182d 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -469,8 +469,8 @@ def _transform_disjunctionDatas( initializer=_setup_spawn, initargs=( dill.dumps(instance), - self._config.solver.name, - self._config.solver.options, + dill.dumps(self._config.solver.__class__), + dill.dumps(self._config.solver.options), self._config.use_primal_bound, ), ) @@ -914,7 +914,7 @@ def get_all_M_values(self, model): # Things we call in subprocesses. These can't be member functions, or # else we'd have to pickle `self`, which is problematic. -def _setup_spawn(model, solver_name, solver_options, use_primal_bound): +def _setup_spawn(model, solver_class, solver_options, use_primal_bound): # When using 'spawn' or 'forkserver', Python starts in a new # environment and executes only this file, so we need to manually # ensure necessary plugins are registered (even if the main process @@ -922,7 +922,7 @@ def _setup_spawn(model, solver_name, solver_options, use_primal_bound): import pyomo.environ _thread_local.model = dill.loads(model) - _thread_local.solver = SolverFactory(solver_name, options=solver_options) + _thread_local.solver = dill.loads(solver_class)(options=dill.loads(solver_options)) _thread_local.config_use_primal_bound = use_primal_bound @@ -931,8 +931,8 @@ def _setup_fork(): # remake the solver instead of using the passed argument. All these # processes are copies of the calling thread so the thread-local # still works. - _thread_local.solver = SolverFactory( - _thread_local.solver.name, options=_thread_local.solver.options + _thread_local.solver = _thread_local.solver.__class__( + options=_thread_local.solver.options ) From 9e93f17f83498dbcfa5b4ba540bdce8a02ff1665 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Jun 2025 14:39:07 -0600 Subject: [PATCH 22/26] Move thread pool setup out of _transform_DisjunctionDatas --- pyomo/gdp/plugins/multiple_bigm.py | 81 +++++++++++++++--------------- 1 file changed, 41 insertions(+), 40 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index f1a8138182d..6c14b16111a 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -444,46 +444,7 @@ def _transform_disjunctionDatas( else os.cpu_count() - 1 ) if threads > 1: - method = ( - self._config.process_start_method - if self._config.process_start_method is not None - else ( - 'spawn' - if os.name == 'nt' - else 'forkserver' if len(threading.enumerate()) > 1 else 'fork' - ) - ) - logger.info( - f"Running {len(jobs)} jobs on {threads} worker " - f"processes with process start method {method}." - ) - if method == 'spawn' or method == 'forkserver': - if not dill_available: - raise GDP_Error( - "Dill is required when spawning processes using " - "methods 'spawn' or 'forkserver', but it could " - "not be imported." - ) - pool = multiprocessing.get_context(method).Pool( - processes=threads, - initializer=_setup_spawn, - initargs=( - dill.dumps(instance), - dill.dumps(self._config.solver.__class__), - dill.dumps(self._config.solver.options), - self._config.use_primal_bound, - ), - ) - elif method == 'fork': - _thread_local.model = instance - _thread_local.solver = self._config.solver - _thread_local.config_use_primal_bound = ( - self._config.use_primal_bound - ) - pool = multiprocessing.get_context('fork').Pool( - processes=threads, initializer=_setup_fork, initargs=() - ) - with pool: + with self._setup_pool(threads, instance, len(jobs)) as pool: results = pool.starmap( func=_calc_M, iterable=[ @@ -838,6 +799,46 @@ def _transform_bound_constraints(self, active_disjuncts, transBlock, Ms): return transformed_constraints + def _setup_pool(self, threads, instance, num_jobs): + method = ( + self._config.process_start_method + if self._config.process_start_method is not None + else ( + 'spawn' + if os.name == 'nt' + else 'forkserver' if len(threading.enumerate()) > 1 else 'fork' + ) + ) + logger.info( + f"Running {num_jobs} jobs on {threads} worker " + f"processes with process start method {method}." + ) + if method == 'spawn' or method == 'forkserver': + if not dill_available: + raise GDP_Error( + "Dill is required when spawning processes using " + "methods 'spawn' or 'forkserver', but it could " + "not be imported." + ) + pool = multiprocessing.get_context(method).Pool( + processes=threads, + initializer=_setup_spawn, + initargs=( + dill.dumps(instance), + dill.dumps(self._config.solver.__class__), + dill.dumps(self._config.solver.options), + self._config.use_primal_bound, + ), + ) + elif method == 'fork': + _thread_local.model = instance + _thread_local.solver = self._config.solver + _thread_local.config_use_primal_bound = self._config.use_primal_bound + pool = multiprocessing.get_context('fork').Pool( + processes=threads, initializer=_setup_fork, initargs=() + ) + return pool + def _add_transformation_block(self, to_block): transBlock, new_block = super()._add_transformation_block(to_block) From c37d30501615532642f163ec757e21bc30ad220f Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 30 Jun 2025 09:47:20 -0600 Subject: [PATCH 23/26] check for a rather implausible error case --- pyomo/gdp/plugins/multiple_bigm.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 6c14b16111a..6cb500ac2b8 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -78,6 +78,7 @@ _thread_local.solver = None _thread_local.model = None _thread_local.config_use_primal_bound = None +_thread_local.in_progress = False def Solver(val): @@ -293,6 +294,11 @@ def __init__(self): self.handlers[Suffix] = self._warn_for_active_suffix def _apply_to(self, instance, **kwds): + # check for the rather implausible error case that + # solver.solve() is a metasolver that indirectly calls this + # transformation again + if _thread_local.in_progress: + raise GDP_Error("gdp.mbigm transformation cannot be called recursively") self.used_args = ComponentMap() with PauseGC(): try: @@ -306,8 +312,10 @@ def _apply_to(self, instance, **kwds): _thread_local.model = None _thread_local.solver = None _thread_local.config_use_primal_bound = None + _thread_local.in_progress = False def _apply_to_impl(self, instance, **kwds): + _thread_local.in_progress = True self._process_arguments(instance, **kwds) if self._config.assume_fixed_vars_permanent: self._bound_visitor.use_fixed_var_values_as_bounds = True From a65467eb5b3d612cc043b48b896ab1bca17f6998 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 30 Jun 2025 09:49:32 -0600 Subject: [PATCH 24/26] terminate subprocesses more gently in an attempt to make codecov work --- pyomo/gdp/plugins/multiple_bigm.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 6cb500ac2b8..faf8370ebd9 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -470,6 +470,8 @@ def _transform_disjunctionDatas( ) in jobs ], ) + pool.close() + pool.join() else: _thread_local.model = instance _thread_local.solver = self._config.solver From 56d12edb0bc78c8da5b04b887c81c4eea11da59b Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 30 Jun 2025 13:32:34 -0600 Subject: [PATCH 25/26] reduce lines of code --- pyomo/gdp/plugins/multiple_bigm.py | 50 ++++++++++-------------------- 1 file changed, 16 insertions(+), 34 deletions(-) diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index faf8370ebd9..663a1c49e39 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -444,6 +444,20 @@ def _transform_disjunctionDatas( ) # (Now exiting the DisjunctionDatas loop) if jobs: + jobs_by_name = [ + ( + constraint.getname(fully_qualified=True), + other_disjunct.getname(fully_qualified=True), + unsuccessful_solve_msg, + is_upper, + ) + for ( + constraint, + other_disjunct, + unsuccessful_solve_msg, + is_upper, + ) in jobs + ] threads = ( self._config.threads if self._config.threads is not None @@ -453,23 +467,7 @@ def _transform_disjunctionDatas( ) if threads > 1: with self._setup_pool(threads, instance, len(jobs)) as pool: - results = pool.starmap( - func=_calc_M, - iterable=[ - ( - constraint.getname(fully_qualified=True), - other_disjunct.getname(fully_qualified=True), - unsuccessful_solve_msg, - is_upper, - ) - for ( - constraint, - other_disjunct, - unsuccessful_solve_msg, - is_upper, - ) in jobs - ], - ) + results = pool.starmap(func=_calc_M, iterable=jobs_by_name) pool.close() pool.join() else: @@ -477,23 +475,7 @@ def _transform_disjunctionDatas( _thread_local.solver = self._config.solver _thread_local.config_use_primal_bound = self._config.use_primal_bound logger.info(f"Running {len(jobs)} jobs single-threaded.") - results = itertools.starmap( - _calc_M, - [ - ( - constraint.getname(fully_qualified=True), - other_disjunct.getname(fully_qualified=True), - unsuccessful_solve_msg, - is_upper, - ) - for ( - constraint, - other_disjunct, - unsuccessful_solve_msg, - is_upper, - ) in jobs - ], - ) + results = itertools.starmap(_calc_M, jobs_by_name) deactivated = set() for (constraint, other_disjunct, _, is_upper), ( M, From e14f0e2046bb3a38b5fc503503dfa65558f30573 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 22 Jul 2025 10:59:46 -0600 Subject: [PATCH 26/26] Address review comments --- pyomo/gdp/plugins/bigm_mixin.py | 2 +- pyomo/gdp/plugins/multiple_bigm.py | 68 +++++++++++++++++++++--------- 2 files changed, 48 insertions(+), 22 deletions(-) diff --git a/pyomo/gdp/plugins/bigm_mixin.py b/pyomo/gdp/plugins/bigm_mixin.py index dc9a1cbf42d..db196d8b72f 100644 --- a/pyomo/gdp/plugins/bigm_mixin.py +++ b/pyomo/gdp/plugins/bigm_mixin.py @@ -17,7 +17,7 @@ import pyomo.contrib.fbbt.interval as interval from pyomo.core import Suffix -logger = logging.getLogger('pyomo.gdp.bigm') +logger = logging.getLogger(__name__) def _convert_M_to_tuple(M, constraint, disjunct=None): diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index 663a1c49e39..fb074128006 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -15,9 +15,10 @@ import multiprocessing import os import threading +import enum from pyomo.common.collections import ComponentMap, ComponentSet -from pyomo.common.config import ConfigDict, ConfigValue +from pyomo.common.config import ConfigDict, ConfigValue, InEnum, PositiveInt from pyomo.common.gc_manager import PauseGC from pyomo.common.modeling import unique_component_name from pyomo.common.dependencies import dill, dill_available @@ -94,11 +95,10 @@ def Solver(val): return val -def ProcessStartMethod(val): - if val in {'spawn', 'fork', 'forkserver', None}: - return val - else: - raise ValueError("Expected one of 'spawn', 'fork', or 'forkserver', or None.") +class ProcessStartMethod(str, enum.Enum): + spawn = 'spawn' + fork = 'fork' + forkserver = 'forkserver' @TransformationFactory.register( @@ -242,11 +242,11 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): 'threads', ConfigValue( default=None, - domain=int, + domain=PositiveInt, description="Number of worker processes to use when estimating M values", doc=""" If not specified, use up to the number of available cores, minus - one. If set to one or less, do not spawn processes, and revert to + one. If set to 1, do not spawn processes, and revert to single-threaded operation. """, ), @@ -255,18 +255,19 @@ class MultipleBigMTransformation(GDP_to_MIP_Transformation, _BigM_MixIn): 'process_start_method', ConfigValue( default=None, - domain=ProcessStartMethod, + domain=InEnum(ProcessStartMethod), description="Start method used for spawning processes during M calculation", doc=""" - Options are 'fork', 'spawn', and 'forkserver', or None. See the Python + Options are the elements of the enum ProcessStartMethod, or equivalently the + strings 'fork', 'spawn', or 'forkserver', or None. See the Python multiprocessing documentation for a full description of each of these. When None is passed, we determine a reasonable default. On POSIX, the default is 'fork', unless we detect that Python has multiple threads at the time the process pool is created, in which case we instead use 'forkserver'. On Windows, the default and only possible option is 'spawn'. Note that if 'spawn' or 'forkserver' are selected, we depend on the `dill` module for - pickling, and model instances must be pickleable using `dill`. This option - is ignored if `threads` is set to one or less. + pickling, and model instances must be pickleable using `dill`. This option is + ignored if `threads` is set to 1. """, ), ) @@ -417,7 +418,8 @@ def _transform_disjunctionDatas( disjunction_setup[t] = (transBlock, algebraic_constraint) = ( self._setup_transform_disjunctionData(t, gdp_tree.root_disjunct(t)) ) - + # Unlike python set(), ComponentSet keeps a stable + # ordering, so we use it for the sake of determinism. active_disjuncts[t] = ComponentSet( disj for disj in t.disjuncts if disj.active ) @@ -532,9 +534,6 @@ def _transform_disjunctionDatas( ) disjunction.deactivate() - # Do the inner work setting up M calculation jobs for one - # disjunction. This mutates the parameters Ms, jobs, - # transBlock._mbm_values, and also self.used_args. def _setup_jobs_for_disjunction( self, disjunction, @@ -545,6 +544,26 @@ def _setup_jobs_for_disjunction( jobs, transBlock, ): + """ + Do the inner work setting up or skipping M calculation jobs for a + disjunction. Mutate the parameters Ms, jobs, transBlock._mbm_values, + and self.used_args. + + Args: + self: self.used_args map from (constraint, disjunct) to M tuple + updated with the keys used from arg_Ms + disjunction: disjunction to set up the jobs for + active_disjuncts: map from disjunctions to ComponentSets of active + disjuncts + arg_Ms: user-provided map from (constraint, disjunct) to M value + or tuple + Ms: working map from (constraint, disjunct) to M tuples to update + jobs: working list of (constraint, other_disjunct, + unsuccessful_solve_msg, is_upper) job tuples to update + transBlock: working transformation block. Update the + transBlock._mbigm_values map from (constraint, disjunct) to + M tuples + """ for disjunct, other_disjunct in itertools.product( active_disjuncts[disjunction], active_disjuncts[disjunction] ): @@ -796,23 +815,30 @@ def _setup_pool(self, threads, instance, num_jobs): self._config.process_start_method if self._config.process_start_method is not None else ( - 'spawn' + ProcessStartMethod.spawn if os.name == 'nt' - else 'forkserver' if len(threading.enumerate()) > 1 else 'fork' + else ( + ProcessStartMethod.forkserver + if len(threading.enumerate()) > 1 + else ProcessStartMethod.fork + ) ) ) logger.info( f"Running {num_jobs} jobs on {threads} worker " f"processes with process start method {method}." ) - if method == 'spawn' or method == 'forkserver': + if ( + method == ProcessStartMethod.spawn + or method == ProcessStartMethod.forkserver + ): if not dill_available: raise GDP_Error( "Dill is required when spawning processes using " "methods 'spawn' or 'forkserver', but it could " "not be imported." ) - pool = multiprocessing.get_context(method).Pool( + pool = multiprocessing.get_context(method.value).Pool( processes=threads, initializer=_setup_spawn, initargs=( @@ -822,7 +848,7 @@ def _setup_pool(self, threads, instance, num_jobs): self._config.use_primal_bound, ), ) - elif method == 'fork': + elif method == ProcessStartMethod.fork: _thread_local.model = instance _thread_local.solver = self._config.solver _thread_local.config_use_primal_bound = self._config.use_primal_bound