1# ___________________________________________________________________________ 2# 3# Pyomo: Python Optimization Modeling Objects 4# Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC 5# Under the terms of Contract DE-NA0003525 with National Technology and 6# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain 7# rights in this software. 8# This software is distributed under the 3-clause BSD License. 9# ___________________________________________________________________________ 10 11"""Big-M Generalized Disjunctive Programming transformation module.""" 12 13import logging 14 15from pyomo.common.collections import ComponentMap, ComponentSet 16from pyomo.common.config import ConfigBlock, ConfigValue 17from pyomo.common.log import is_debug_set 18from pyomo.common.modeling import unique_component_name 19from pyomo.common.deprecation import deprecated, deprecation_warning 20from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr 21from pyomo.core import ( 22 Block, BooleanVar, Connector, Constraint, Param, Set, SetOf, Suffix, Var, 23 Expression, SortComponents, TraversalStrategy, value, 24 RangeSet, NonNegativeIntegers, LogicalConstraint, ) 25from pyomo.core.base.external import ExternalFunction 26from pyomo.core.base import Transformation, TransformationFactory, Reference 27import pyomo.core.expr.current as EXPR 28from pyomo.gdp import Disjunct, Disjunction, GDP_Error 29from pyomo.gdp.util import ( _warn_for_active_logical_constraint, target_list, 30 is_child_of, get_src_disjunction, 31 get_src_constraint, get_transformed_constraints, 32 _get_constraint_transBlock, get_src_disjunct, 33 _warn_for_active_disjunction, 34 _warn_for_active_disjunct, preprocess_targets) 35from pyomo.network import Port 36from pyomo.repn import generate_standard_repn 37from functools import wraps 38from weakref import ref as weakref_ref 39 40logger = logging.getLogger('pyomo.gdp.bigm') 41 42NAME_BUFFER = {} 43 44def _to_dict(val): 45 if isinstance(val, (dict, ComponentMap)): 46 return val 47 return {None: val} 48 49 50@TransformationFactory.register('gdp.bigm', doc="Relax disjunctive model using " 51 "big-M terms.") 52class BigM_Transformation(Transformation): 53 """Relax disjunctive model using big-M terms. 54 55 Relaxes a disjunctive model into an algebraic model by adding Big-M 56 terms to all disjunctive constraints. 57 58 This transformation accepts the following keyword arguments: 59 bigM: A user-specified value (or dict) of M values to use (see below) 60 targets: the targets to transform [default: the instance] 61 62 M values are determined as follows: 63 1) if the constraint appears in the bigM argument dict 64 2) if the constraint parent_component appears in the bigM 65 argument dict 66 3) if any block which is an ancestor to the constraint appears in 67 the bigM argument dict 68 3) if 'None' is in the bigM argument dict 69 4) if the constraint or the constraint parent_component appear in 70 a BigM Suffix attached to any parent_block() beginning with the 71 constraint's parent_block and moving up to the root model. 72 5) if None appears in a BigM Suffix attached to any 73 parent_block() between the constraint and the root model. 74 6) if the constraint is linear, estimate M using the variable bounds 75 76 M values may be a single value or a 2-tuple specifying the M for the 77 lower bound and the upper bound of the constraint body. 78 79 Specifying "bigM=N" is automatically mapped to "bigM={None: N}". 80 81 The transformation will create a new Block with a unique 82 name beginning "_pyomo_gdp_bigm_reformulation". That Block will 83 contain an indexed Block named "relaxedDisjuncts", which will hold 84 the relaxed disjuncts. This block is indexed by an integer 85 indicating the order in which the disjuncts were relaxed. 86 Each block has a dictionary "_constraintMap": 87 88 'srcConstraints': ComponentMap(<transformed constraint>: 89 <src constraint>) 90 'transformedConstraints': ComponentMap(<src constraint>: 91 <transformed constraint>) 92 93 All transformed Disjuncts will have a pointer to the block their transformed 94 constraints are on, and all transformed Disjunctions will have a 95 pointer to the corresponding OR or XOR constraint. 96 97 """ 98 99 CONFIG = ConfigBlock("gdp.bigm") 100 CONFIG.declare('targets', ConfigValue( 101 default=None, 102 domain=target_list, 103 description="target or list of targets that will be relaxed", 104 doc=""" 105 106 This specifies the list of components to relax. If None (default), the 107 entire model is transformed. Note that if the transformation is done out 108 of place, the list of targets should be attached to the model before it 109 is cloned, and the list will specify the targets on the cloned 110 instance.""" 111 )) 112 CONFIG.declare('bigM', ConfigValue( 113 default=None, 114 domain=_to_dict, 115 description="Big-M value used for constraint relaxation", 116 doc=""" 117 118 A user-specified value, dict, or ComponentMap of M values that override 119 M-values found through model Suffixes or that would otherwise be 120 calculated using variable domains.""" 121 )) 122 CONFIG.declare('assume_fixed_vars_permanent', ConfigValue( 123 default=False, 124 domain=bool, 125 description="Boolean indicating whether or not to transform so that the " 126 "the transformed model will still be valid when fixed Vars are unfixed.", 127 doc=""" 128 This is only relevant when the transformation will be estimating values 129 for M. If True, the transformation will calculate M values assuming that 130 fixed variables will always be fixed to their current values. This means 131 that if a fixed variable is unfixed after transformation, the 132 transformed model is potentially no longer valid. By default, the 133 transformation will assume fixed variables could be unfixed in the 134 future and will use their bounds to calculate the M value rather than 135 their value. Note that this could make for a weaker LP relaxation 136 while the variables remain fixed. 137 """ 138 )) 139 140 def __init__(self): 141 """Initialize transformation object.""" 142 super(BigM_Transformation, self).__init__() 143 self.handlers = { 144 Constraint: self._transform_constraint, 145 Var: False, # Note that if a Var appears on a Disjunct, we 146 # still treat its bounds as global. If the 147 # intent is for its bounds to be on the 148 # disjunct, it should be declared with no bounds 149 # and the bounds should be set in constraints on 150 # the Disjunct. 151 BooleanVar: False, 152 Connector: False, 153 Expression: False, 154 Suffix: False, 155 Param: False, 156 Set: False, 157 SetOf: False, 158 RangeSet: False, 159 Disjunction: self._warn_for_active_disjunction, 160 Disjunct: self._warn_for_active_disjunct, 161 Block: self._transform_block_on_disjunct, 162 LogicalConstraint: self._warn_for_active_logical_statement, 163 ExternalFunction: False, 164 Port: False, # not Arcs, because those are deactivated after 165 # the network.expand_arcs transformation 166 } 167 self._generate_debug_messages = False 168 169 def _get_bigm_suffix_list(self, block, stopping_block=None): 170 # Note that you can only specify suffixes on BlockData objects or 171 # ScalarBlocks. Though it is possible at this point to stick them 172 # on whatever components you want, we won't pick them up. 173 suffix_list = [] 174 175 # go searching above block in the tree, stop when we hit stopping_block 176 # (This is so that we can search on each Disjunct once, but get any 177 # information between a constraint and its Disjunct while transforming 178 # the constraint). 179 while block is not stopping_block: 180 bigm = block.component('BigM') 181 if type(bigm) is Suffix: 182 suffix_list.append(bigm) 183 block = block.parent_block() 184 185 return suffix_list 186 187 def _get_bigm_arg_list(self, bigm_args, block): 188 # Gather what we know about blocks from args exactly once. We'll still 189 # check for constraints in the moment, but if that fails, we've 190 # preprocessed the time-consuming part of traversing up the tree. 191 arg_list = [] 192 if bigm_args is None: 193 return arg_list 194 while block is not None: 195 if block in bigm_args: 196 arg_list.append({block: bigm_args[block]}) 197 block = block.parent_block() 198 return arg_list 199 200 def _apply_to(self, instance, **kwds): 201 assert not NAME_BUFFER 202 self._generate_debug_messages = is_debug_set(logger) 203 self.used_args = ComponentMap() # If everything was sure to go well, 204 # this could be a dictionary. But if 205 # someone messes up and gives us a Var 206 # as a key in bigMargs, I need the error 207 # not to be when I try to put it into 208 # this map! 209 try: 210 self._apply_to_impl(instance, **kwds) 211 finally: 212 # Clear the global name buffer now that we are done 213 NAME_BUFFER.clear() 214 # same for our bookkeeping about what we used from bigM arg dict 215 self.used_args.clear() 216 217 def _apply_to_impl(self, instance, **kwds): 218 if not instance.ctype in (Block, Disjunct): 219 raise GDP_Error("Transformation called on %s of type %s. 'instance' " 220 "must be a ConcreteModel, Block, or Disjunct (in " 221 "the case of nested disjunctions)." % 222 (instance.name, instance.ctype)) 223 224 config = self.CONFIG(kwds.pop('options', {})) 225 226 # We will let args override suffixes and estimate as a last 227 # resort. More specific args/suffixes override ones anywhere in 228 # the tree. Suffixes lower down in the tree override ones higher 229 # up. 230 if 'default_bigM' in kwds: 231 deprecation_warning("the 'default_bigM=' argument has been " 232 "replaced by 'bigM='", version='5.4') 233 config.bigM = kwds.pop('default_bigM') 234 235 config.set_value(kwds) 236 bigM = config.bigM 237 self.assume_fixed_vars_permanent = config.assume_fixed_vars_permanent 238 239 targets = config.targets 240 if targets is None: 241 targets = (instance, ) 242 else: 243 # we need to preprocess targets to make sure that if there are any 244 # disjunctions in targets that their disjuncts appear before them in 245 # the list. 246 targets = preprocess_targets(targets) 247 248 # We need to check that all the targets are in fact on 249 # instance. As we do this, we will use the set below to cache components 250 # we know to be in the tree rooted at instance. 251 knownBlocks = {} 252 for t in targets: 253 # check that t is in fact a child of instance 254 if not is_child_of(parent=instance, child=t, 255 knownBlocks=knownBlocks): 256 raise GDP_Error( 257 "Target '%s' is not a component on instance '%s'!" 258 % (t.name, instance.name)) 259 elif t.ctype is Disjunction: 260 if t.is_indexed(): 261 self._transform_disjunction(t, bigM) 262 else: 263 self._transform_disjunctionData( t, bigM, t.index()) 264 elif t.ctype in (Block, Disjunct): 265 if t.is_indexed(): 266 self._transform_block(t, bigM) 267 else: 268 self._transform_blockData(t, bigM) 269 else: 270 raise GDP_Error( 271 "Target '%s' was not a Block, Disjunct, or Disjunction. " 272 "It was of type %s and can't be transformed." 273 % (t.name, type(t))) 274 275 # issue warnings about anything that was in the bigM args dict that we 276 # didn't use 277 if bigM is not None: 278 unused_args = ComponentSet(bigM.keys()) - \ 279 ComponentSet(self.used_args.keys()) 280 if len(unused_args) > 0: 281 warning_msg = ("Unused arguments in the bigM map! " 282 "These arguments were not used by the " 283 "transformation:\n") 284 for component in unused_args: 285 if hasattr(component, 'name'): 286 warning_msg += "\t%s\n" % component.name 287 else: 288 warning_msg += "\t%s\n" % component 289 logger.warning(warning_msg) 290 291 def _add_transformation_block(self, instance): 292 # make a transformation block on instance to put transformed disjuncts 293 # on 294 transBlockName = unique_component_name( 295 instance, 296 '_pyomo_gdp_bigm_reformulation') 297 transBlock = Block() 298 instance.add_component(transBlockName, transBlock) 299 transBlock.relaxedDisjuncts = Block(NonNegativeIntegers) 300 transBlock.lbub = Set(initialize=['lb', 'ub']) 301 302 return transBlock 303 304 def _transform_block(self, obj, bigM): 305 for i in sorted(obj.keys()): 306 self._transform_blockData(obj[i], bigM) 307 308 def _transform_blockData(self, obj, bigM): 309 # Transform every (active) disjunction in the block 310 for disjunction in obj.component_objects( 311 Disjunction, 312 active=True, 313 sort=SortComponents.deterministic, 314 descend_into=(Block, Disjunct), 315 descent_order=TraversalStrategy.PostfixDFS): 316 self._transform_disjunction(disjunction, bigM) 317 318 def _add_xor_constraint(self, disjunction, transBlock): 319 # Put the disjunction constraint on the transformation block and 320 # determine whether it is an OR or XOR constraint. 321 322 # We never do this for just a DisjunctionData because we need to know 323 # about the index set of its parent component (so that we can make the 324 # index of this constraint match). So if we called this on a 325 # DisjunctionData, we did something wrong. 326 assert isinstance(disjunction, Disjunction) 327 328 # first check if the constraint already exists 329 if disjunction._algebraic_constraint is not None: 330 return disjunction._algebraic_constraint() 331 332 # add the XOR (or OR) constraints to parent block (with unique name) 333 # It's indexed if this is an IndexedDisjunction, not otherwise 334 orC = Constraint(disjunction.index_set()) if \ 335 disjunction.is_indexed() else Constraint() 336 # The name used to indicate if there were OR or XOR disjunctions, 337 # however now that Disjunctions are allowed to mix the state we 338 # can no longer make that distinction in the name. 339 # nm = '_xor' if xor else '_or' 340 nm = '_xor' 341 orCname = unique_component_name( transBlock, disjunction.getname( 342 fully_qualified=True, name_buffer=NAME_BUFFER) + nm) 343 transBlock.add_component(orCname, orC) 344 disjunction._algebraic_constraint = weakref_ref(orC) 345 346 return orC 347 348 def _transform_disjunction(self, obj, bigM): 349 if not obj.active: 350 return 351 352 # if this is an IndexedDisjunction we have seen in a prior call to the 353 # transformation, we already have a transformation block for it. We'll 354 # use that. 355 if obj._algebraic_constraint is not None: 356 transBlock = obj._algebraic_constraint().parent_block() 357 else: 358 transBlock = self._add_transformation_block(obj.parent_block()) 359 360 # relax each of the disjunctionDatas 361 for i in sorted(obj.keys()): 362 self._transform_disjunctionData(obj[i], bigM, i, transBlock) 363 364 # deactivate so the writers don't scream 365 obj.deactivate() 366 367 def _transform_disjunctionData(self, obj, bigM, index, transBlock=None): 368 if not obj.active: 369 return # Do not process a deactivated disjunction 370 # We won't have these arguments if this got called straight from 371 # targets. But else, we created them earlier, and have just been passing 372 # them through. 373 if transBlock is None: 374 # It's possible that we have already created a transformation block 375 # for another disjunctionData from this same container. If that's 376 # the case, let's use the same transformation block. (Else it will 377 # be really confusing that the XOR constraint goes to that old block 378 # but we create a new one here.) 379 if obj.parent_component()._algebraic_constraint is not None: 380 transBlock = obj.parent_component()._algebraic_constraint().\ 381 parent_block() 382 else: 383 transBlock = self._add_transformation_block(obj.parent_block()) 384 # create or fetch the xor constraint 385 xorConstraint = self._add_xor_constraint(obj.parent_component(), 386 transBlock) 387 388 xor = obj.xor 389 or_expr = 0 390 # Just because it's unlikely this is what someone meant to do... 391 if len(obj.disjuncts) == 0: 392 raise GDP_Error("Disjunction '%s' is empty. This is " 393 "likely indicative of a modeling error." % 394 obj.getname(fully_qualified=True, 395 name_buffer=NAME_BUFFER)) 396 for disjunct in obj.disjuncts: 397 or_expr += disjunct.binary_indicator_var 398 # make suffix list. (We don't need it until we are 399 # transforming constraints, but it gets created at the 400 # disjunct level, so more efficient to make it here and 401 # pass it down.) 402 suffix_list = self._get_bigm_suffix_list(disjunct) 403 arg_list = self._get_bigm_arg_list(bigM, disjunct) 404 # relax the disjunct 405 self._transform_disjunct(disjunct, transBlock, bigM, arg_list, 406 suffix_list) 407 408 # add or (or xor) constraint 409 if xor: 410 xorConstraint[index] = or_expr == 1 411 else: 412 xorConstraint[index] = or_expr >= 1 413 # Mark the DisjunctionData as transformed by mapping it to its XOR 414 # constraint. 415 obj._algebraic_constraint = weakref_ref(xorConstraint[index]) 416 417 # and deactivate for the writers 418 obj.deactivate() 419 420 def _transform_disjunct(self, obj, transBlock, bigM, arg_list, suffix_list): 421 # deactivated -> either we've already transformed or user deactivated 422 if not obj.active: 423 if obj.indicator_var.is_fixed(): 424 if not value(obj.indicator_var): 425 # The user cleanly deactivated the disjunct: there 426 # is nothing for us to do here. 427 return 428 else: 429 raise GDP_Error( 430 "The disjunct '%s' is deactivated, but the " 431 "indicator_var is fixed to %s. This makes no sense." 432 % ( obj.name, value(obj.indicator_var) )) 433 if obj._transformation_block is None: 434 raise GDP_Error( 435 "The disjunct '%s' is deactivated, but the " 436 "indicator_var is not fixed and the disjunct does not " 437 "appear to have been relaxed. This makes no sense. " 438 "(If the intent is to deactivate the disjunct, fix its " 439 "indicator_var to False.)" 440 % ( obj.name, )) 441 442 if obj._transformation_block is not None: 443 # we've transformed it, which means this is the second time it's 444 # appearing in a Disjunction 445 raise GDP_Error( 446 "The disjunct '%s' has been transformed, but a disjunction " 447 "it appears in has not. Putting the same disjunct in " 448 "multiple disjunctions is not supported." % obj.name) 449 450 # add reference to original disjunct on transformation block 451 relaxedDisjuncts = transBlock.relaxedDisjuncts 452 relaxationBlock = relaxedDisjuncts[len(relaxedDisjuncts)] 453 # we will keep a map of constraints (hashable, ha!) to a tuple to 454 # indicate what their M value is and where it came from, of the form: 455 # ((lower_value, lower_source, lower_key), (upper_value, upper_source, 456 # upper_key)), where the first tuple is the information for the lower M, 457 # the second tuple is the info for the upper M, source is the Suffix or 458 # argument dictionary and None if the value was calculated, and key is 459 # the key in the Suffix or argument dictionary, and None if it was 460 # calculated. (Note that it is possible the lower or upper is 461 # user-specified and the other is not, hence the need to store 462 # information for both.) 463 relaxationBlock.bigm_src = {} 464 relaxationBlock.localVarReferences = Block() 465 obj._transformation_block = weakref_ref(relaxationBlock) 466 relaxationBlock._srcDisjunct = weakref_ref(obj) 467 468 # This is crazy, but if the disjunction has been previously 469 # relaxed, the disjunct *could* be deactivated. This is a big 470 # deal for Hull, as it uses the component_objects / 471 # component_data_objects generators. For BigM, that is OK, 472 # because we never use those generators with active=True. I am 473 # only noting it here for the future when someone (me?) is 474 # comparing the two relaxations. 475 # 476 # Transform each component within this disjunct 477 self._transform_block_components(obj, obj, bigM, arg_list, suffix_list) 478 479 # deactivate disjunct to keep the writers happy 480 obj._deactivate_without_fixing_indicator() 481 482 def _transform_block_components(self, block, disjunct, bigM, arg_list, 483 suffix_list): 484 # We find any transformed disjunctions that might be here because we 485 # need to move their transformation blocks up onto the parent block 486 # before we transform anything else on this block. Note that we do this 487 # before we create references to local variables because we do not want 488 # duplicate references to indicator variables and local variables on 489 # nested disjuncts. 490 disjunctBlock = disjunct._transformation_block() 491 destinationBlock = disjunctBlock.parent_block() 492 for obj in block.component_data_objects( 493 Disjunction, 494 sort=SortComponents.deterministic, 495 descend_into=(Block)): 496 if obj.algebraic_constraint is None: 497 # This could be bad if it's active since that means its 498 # untransformed, but we'll wait to yell until the next loop 499 continue 500 # get this disjunction's relaxation block. 501 transBlock = obj.algebraic_constraint().parent_block() 502 503 # move transBlock up to parent component 504 self._transfer_transBlock_data(transBlock, destinationBlock) 505 # we leave the transformation block because it still has the XOR 506 # constraints, which we want to be on the parent disjunct. 507 508 # Find all the variables declared here (including the indicator_var) and 509 # add a reference on the transformation block so these will be 510 # accessible when the Disjunct is deactivated. We don't descend into 511 # Disjuncts because we just moved the references to their local 512 # variables up in the previous loop. 513 varRefBlock = disjunctBlock.localVarReferences 514 for v in block.component_objects(Var, descend_into=Block, active=None): 515 varRefBlock.add_component(unique_component_name( 516 varRefBlock, v.getname(fully_qualified=True, 517 name_buffer=NAME_BUFFER)), Reference(v)) 518 519 # Now look through the component map of block and transform everything 520 # we have a handler for. Yell if we don't know how to handle it. (Note 521 # that because we only iterate through active components, this means 522 # non-ActiveComponent types cannot have handlers.) 523 for obj in block.component_objects(active=True, descend_into=False): 524 handler = self.handlers.get(obj.ctype, None) 525 if not handler: 526 if handler is None: 527 raise GDP_Error( 528 "No BigM transformation handler registered " 529 "for modeling components of type %s. If your " 530 "disjuncts contain non-GDP Pyomo components that " 531 "require transformation, please transform them first." 532 % obj.ctype) 533 continue 534 # obj is what we are transforming, we pass disjunct 535 # through so that we will have access to the indicator 536 # variables down the line. 537 handler(obj, disjunct, bigM, arg_list, suffix_list) 538 539 def _transfer_transBlock_data(self, fromBlock, toBlock): 540 # We know that we have a list of transformed disjuncts on both. We need 541 # to move those over. We know the XOR constraints are on the block, and 542 # we need to leave those on the disjunct. 543 disjunctList = toBlock.relaxedDisjuncts 544 to_delete = [] 545 for idx, disjunctBlock in fromBlock.relaxedDisjuncts.items(): 546 newblock = disjunctList[len(disjunctList)] 547 newblock.transfer_attributes_from(disjunctBlock) 548 549 # update the mappings 550 original = disjunctBlock._srcDisjunct() 551 original._transformation_block = weakref_ref(newblock) 552 newblock._srcDisjunct = weakref_ref(original) 553 554 # save index of what we just moved so that we can delete it 555 to_delete.append(idx) 556 557 # delete everything we moved. 558 for idx in to_delete: 559 del fromBlock.relaxedDisjuncts[idx] 560 561 # Note that we could handle other components here if we ever needed 562 # to, but we control what is on the transformation block and 563 # currently everything is on the blocks that we just moved... 564 565 def _warn_for_active_disjunction(self, disjunction, disjunct, bigMargs, 566 arg_list, suffix_list): 567 _warn_for_active_disjunction(disjunction, disjunct, NAME_BUFFER) 568 569 def _warn_for_active_disjunct(self, innerdisjunct, outerdisjunct, bigMargs, 570 arg_list, suffix_list): 571 _warn_for_active_disjunct(innerdisjunct, outerdisjunct, NAME_BUFFER) 572 573 def _warn_for_active_logical_statement( 574 self, logical_statment, disjunct, infodict, bigMargs, suffix_list): 575 _warn_for_active_logical_constraint(logical_statment, disjunct, 576 NAME_BUFFER) 577 578 def _transform_block_on_disjunct(self, block, disjunct, bigMargs, arg_list, 579 suffix_list): 580 # We look through everything on the component map of the block 581 # and transform it just as we would if it was on the disjunct 582 # directly. (We are passing the disjunct through so that when 583 # we find constraints, _xform_constraint will have access to 584 # the correct indicator variable.) 585 for i in sorted(block.keys()): 586 self._transform_block_components( block[i], disjunct, bigMargs, 587 arg_list, suffix_list) 588 589 def _get_constraint_map_dict(self, transBlock): 590 if not hasattr(transBlock, "_constraintMap"): 591 transBlock._constraintMap = { 592 'srcConstraints': ComponentMap(), 593 'transformedConstraints': ComponentMap()} 594 return transBlock._constraintMap 595 596 def _convert_M_to_tuple(self, M, constraint_name): 597 if not isinstance(M, (tuple, list)): 598 if M is None: 599 M = (None, None) 600 else: 601 try: 602 M = (-M, M) 603 except: 604 logger.error("Error converting scalar M-value %s " 605 "to (-M,M). Is %s not a numeric type?" 606 % (M, type(M))) 607 raise 608 if len(M) != 2: 609 raise GDP_Error("Big-M %s for constraint %s is not of " 610 "length two. " 611 "Expected either a single value or " 612 "tuple or list of length two for M." 613 % (str(M), constraint_name)) 614 615 return M 616 617 def _transform_constraint(self, obj, disjunct, bigMargs, arg_list, 618 disjunct_suffix_list): 619 # add constraint to the transformation block, we'll transform it there. 620 transBlock = disjunct._transformation_block() 621 bigm_src = transBlock.bigm_src 622 constraintMap = self._get_constraint_map_dict(transBlock) 623 624 disjunctionRelaxationBlock = transBlock.parent_block() 625 # Though rare, it is possible to get naming conflicts here 626 # since constraints from all blocks are getting moved onto the 627 # same block. So we get a unique name 628 cons_name = obj.getname(fully_qualified=True, name_buffer=NAME_BUFFER) 629 name = unique_component_name(transBlock, cons_name) 630 631 if obj.is_indexed(): 632 newConstraint = Constraint(obj.index_set(), 633 disjunctionRelaxationBlock.lbub) 634 # we map the container of the original to the container of the 635 # transformed constraint. Don't do this if obj is a ScalarConstraint 636 # because we will treat that like a _ConstraintData and map to a 637 # list of transformed _ConstraintDatas 638 constraintMap['transformedConstraints'][obj] = newConstraint 639 else: 640 newConstraint = Constraint(disjunctionRelaxationBlock.lbub) 641 transBlock.add_component(name, newConstraint) 642 # add mapping of transformed constraint to original constraint 643 constraintMap['srcConstraints'][newConstraint] = obj 644 645 for i in sorted(obj.keys()): 646 c = obj[i] 647 if not c.active: 648 continue 649 650 lower = (None, None, None) 651 upper = (None, None, None) 652 653 # first, we see if an M value was specified in the arguments. 654 # (This returns None if not) 655 lower, upper = self._get_M_from_args(c, bigMargs, arg_list, lower, 656 upper) 657 M = (lower[0], upper[0]) 658 659 if self._generate_debug_messages: 660 _name = obj.getname( 661 fully_qualified=True, name_buffer=NAME_BUFFER) 662 logger.debug("GDP(BigM): The value for M for constraint '%s' " 663 "from the BigM argument is %s." % (cons_name, 664 str(M))) 665 666 # if we didn't get something we need from args, try suffixes: 667 if (M[0] is None and c.lower is not None) or \ 668 (M[1] is None and c.upper is not None): 669 # first get anything parent to c but below disjunct 670 suffix_list = self._get_bigm_suffix_list(c.parent_block(), 671 stopping_block=disjunct) 672 # prepend that to what we already collected for the disjunct. 673 suffix_list.extend(disjunct_suffix_list) 674 lower, upper = self._update_M_from_suffixes(c, suffix_list, 675 lower, upper) 676 M = (lower[0], upper[0]) 677 678 if self._generate_debug_messages: 679 _name = obj.getname( 680 fully_qualified=True, name_buffer=NAME_BUFFER) 681 logger.debug("GDP(BigM): The value for M for constraint '%s' " 682 "after checking suffixes is %s." % (cons_name, 683 str(M))) 684 685 if c.lower is not None and M[0] is None: 686 M = (self._estimate_M(c.body, name)[0] - c.lower, M[1]) 687 lower = (M[0], None, None) 688 if c.upper is not None and M[1] is None: 689 M = (M[0], self._estimate_M(c.body, name)[1] - c.upper) 690 upper = (M[1], None, None) 691 692 if self._generate_debug_messages: 693 _name = obj.getname( 694 fully_qualified=True, name_buffer=NAME_BUFFER) 695 logger.debug("GDP(BigM): The value for M for constraint '%s' " 696 "after estimating (if needed) is %s." % 697 (cons_name, str(M))) 698 699 # save the source information 700 bigm_src[c] = (lower, upper) 701 702 # Handle indices for both ScalarConstraint and IndexedConstraint 703 if i.__class__ is tuple: 704 i_lb = i + ('lb',) 705 i_ub = i + ('ub',) 706 elif obj.is_indexed(): 707 i_lb = (i, 'lb',) 708 i_ub = (i, 'ub',) 709 else: 710 i_lb = 'lb' 711 i_ub = 'ub' 712 713 if c.lower is not None: 714 if M[0] is None: 715 raise GDP_Error("Cannot relax disjunctive constraint '%s' " 716 "because M is not defined." % name) 717 M_expr = M[0] * (1 - disjunct.binary_indicator_var) 718 newConstraint.add(i_lb, c.lower <= c. body - M_expr) 719 constraintMap[ 720 'transformedConstraints'][c] = [newConstraint[i_lb]] 721 constraintMap['srcConstraints'][newConstraint[i_lb]] = c 722 if c.upper is not None: 723 if M[1] is None: 724 raise GDP_Error("Cannot relax disjunctive constraint '%s' " 725 "because M is not defined." % name) 726 M_expr = M[1] * (1 - disjunct.binary_indicator_var) 727 newConstraint.add(i_ub, c.body - M_expr <= c.upper) 728 transformed = constraintMap['transformedConstraints'].get(c) 729 if transformed is not None: 730 constraintMap['transformedConstraints'][ 731 c].append(newConstraint[i_ub]) 732 else: 733 constraintMap[ 734 'transformedConstraints'][c] = [newConstraint[i_ub]] 735 constraintMap['srcConstraints'][newConstraint[i_ub]] = c 736 737 # deactivate because we relaxed 738 c.deactivate() 739 740 def _process_M_value(self, m, lower, upper, need_lower, need_upper, src, 741 key, constraint_name, from_args=False): 742 m = self._convert_M_to_tuple(m, constraint_name) 743 if need_lower and m[0] is not None: 744 if from_args: 745 self.used_args[key] = m 746 lower = (m[0], src, key) 747 need_lower = False 748 if need_upper and m[1] is not None: 749 if from_args: 750 self.used_args[key] = m 751 upper = (m[1], src, key) 752 need_upper = False 753 return lower, upper, need_lower, need_upper 754 755 def _get_M_from_args(self, constraint, bigMargs, arg_list, lower, upper): 756 # check args: we first look in the keys for constraint and 757 # constraintdata. In the absence of those, we traverse up the blocks, 758 # and as a last resort check for a value for None 759 if bigMargs is None: 760 return (lower, upper) 761 762 # since we check for args first, we know lower[0] and upper[0] are both 763 # None 764 need_lower = constraint.lower is not None 765 need_upper = constraint.upper is not None 766 constraint_name = constraint.getname(fully_qualified=True, 767 name_buffer=NAME_BUFFER) 768 769 # check for the constraint itself and its container 770 parent = constraint.parent_component() 771 if constraint in bigMargs: 772 m = bigMargs[constraint] 773 (lower, upper, 774 need_lower, need_upper) = self._process_M_value(m, lower, upper, 775 need_lower, 776 need_upper, 777 bigMargs, 778 constraint, 779 constraint_name, 780 from_args=True) 781 if not need_lower and not need_upper: 782 return lower, upper 783 elif parent in bigMargs: 784 m = bigMargs[parent] 785 (lower, upper, 786 need_lower, need_upper) = self._process_M_value(m, lower, upper, 787 need_lower, 788 need_upper, 789 bigMargs, parent, 790 constraint_name, 791 from_args=True) 792 if not need_lower and not need_upper: 793 return lower, upper 794 795 # use the precomputed traversal up the blocks 796 for arg in arg_list: 797 for block, val in arg.items(): 798 (lower, upper, 799 need_lower, need_upper) = self._process_M_value(val, lower, 800 upper, 801 need_lower, 802 need_upper, 803 bigMargs, 804 block, 805 constraint_name, 806 from_args=True) 807 if not need_lower and not need_upper: 808 return lower, upper 809 810 # last check for value for None! 811 if None in bigMargs: 812 m = bigMargs[None] 813 (lower, upper, 814 need_lower, need_upper) = self._process_M_value(m, lower, upper, 815 need_lower, 816 need_upper, 817 bigMargs, None, 818 constraint_name, 819 from_args=True) 820 if not need_lower and not need_upper: 821 return lower, upper 822 823 return lower, upper 824 825 def _update_M_from_suffixes(self, constraint, suffix_list, lower, upper): 826 # It's possible we found half the answer in args, but we are still 827 # looking for half the answer. 828 need_lower = constraint.lower is not None and lower[0] is None 829 need_upper = constraint.upper is not None and upper[0] is None 830 constraint_name = constraint.getname(fully_qualified=True, 831 name_buffer=NAME_BUFFER) 832 M = None 833 # first we check if the constraint or its parent is a key in any of the 834 # suffix lists 835 for bigm in suffix_list: 836 if constraint in bigm: 837 M = bigm[constraint] 838 (lower, upper, 839 need_lower, need_upper) = self._process_M_value(M, lower, 840 upper, 841 need_lower, 842 need_upper, 843 bigm, 844 constraint, 845 constraint_name) 846 if not need_lower and not need_upper: 847 return lower, upper 848 849 # if c is indexed, check for the parent component 850 if constraint.parent_component() in bigm: 851 parent = constraint.parent_component() 852 M = bigm[parent] 853 (lower, upper, 854 need_lower, need_upper) = self._process_M_value(M, lower, 855 upper, 856 need_lower, 857 need_upper, 858 bigm, parent, 859 constraint_name) 860 if not need_lower and not need_upper: 861 return lower, upper 862 863 # if we didn't get an M that way, traverse upwards through the blocks 864 # and see if None has a value on any of them. 865 if M is None: 866 for bigm in suffix_list: 867 if None in bigm: 868 M = bigm[None] 869 (lower, upper, 870 need_lower, 871 need_upper) = self._process_M_value(M, lower, upper, 872 need_lower, need_upper, 873 bigm, None, 874 constraint_name) 875 if not need_lower and not need_upper: 876 return lower, upper 877 return lower, upper 878 879 def _estimate_M(self, expr, name): 880 # If there are fixed variables here, unfix them for this calculation, 881 # and we'll restore them at the end. 882 fixed_vars = ComponentMap() 883 if not self.assume_fixed_vars_permanent: 884 for v in EXPR.identify_variables(expr, include_fixed=True): 885 if v.fixed: 886 fixed_vars[v] = value(v) 887 v.fixed = False 888 889 expr_lb, expr_ub = compute_bounds_on_expr(expr) 890 if expr_lb is None or expr_ub is None: 891 raise GDP_Error("Cannot estimate M for unbounded " 892 "expressions.\n\t(found while processing " 893 "constraint '%s'). Please specify a value of M " 894 "or ensure all variables that appear in the " 895 "constraint are bounded." % name) 896 else: 897 M = (expr_lb, expr_ub) 898 899 # clean up if we unfixed things (fixed_vars is empty if we were assuming 900 # fixed vars are fixed for life) 901 for v, val in fixed_vars.items(): 902 v.fix(val) 903 904 return tuple(M) 905 906 # These are all functions to retrieve transformed components from 907 # original ones and vice versa. 908 909 @wraps(get_src_disjunct) 910 def get_src_disjunct(self, transBlock): 911 return get_src_disjunct(transBlock) 912 913 @wraps(get_src_disjunction) 914 def get_src_disjunction(self, xor_constraint): 915 return get_src_disjunction(xor_constraint) 916 917 @wraps(get_src_constraint) 918 def get_src_constraint(self, transformedConstraint): 919 return get_src_constraint(transformedConstraint) 920 921 @wraps(get_transformed_constraints) 922 def get_transformed_constraints(self, srcConstraint): 923 return get_transformed_constraints(srcConstraint) 924 925 @deprecated("The get_m_value_src function is deprecated. Use " 926 "the get_M_value_src function is you need source " 927 "information or the get_M_value function if you " 928 "only need values.", version='5.7.1') 929 def get_m_value_src(self, constraint): 930 transBlock = _get_constraint_transBlock(constraint) 931 ((lower_val, lower_source, lower_key), 932 (upper_val, upper_source, upper_key)) = transBlock.bigm_src[constraint] 933 934 if constraint.lower is not None and constraint.upper is not None and \ 935 (not lower_source is upper_source or not lower_key is upper_key): 936 raise GDP_Error("This is why this method is deprecated: The lower " 937 "and upper M values for constraint %s came from " 938 "different sources, please use the get_M_value_src " 939 "method." % constraint.name) 940 # if source and key are equal for the two, this is representable in the 941 # old format. 942 if constraint.lower is not None and lower_source is not None: 943 return (lower_source, lower_key) 944 if constraint.upper is not None and upper_source is not None: 945 return (upper_source, upper_key) 946 # else it was calculated: 947 return (lower_val, upper_val) 948 949 def get_M_value_src(self, constraint): 950 """Return a tuple indicating how the M value used to transform 951 constraint was specified. (In particular, this can be used to 952 verify which BigM Suffixes were actually necessary to the 953 transformation.) 954 955 Return is of the form: ((lower_M_val, lower_M_source, lower_M_key), 956 (upper_M_val, upper_M_source, upper_M_key)) 957 958 If the constraint does not have a lower bound (or an upper bound), 959 the first (second) element will be (None, None, None). Note that if 960 a constraint is of the form a <= expr <= b or is an equality constraint, 961 it is not necessarily true that the source of lower_M and upper_M 962 are the same. 963 964 If the M value came from an arg, source is the dictionary itself and 965 key is the key in that dictionary which gave us the M value. 966 967 If the M value came from a Suffix, source is the BigM suffix used and 968 key is the key in that Suffix. 969 970 If the transformation calculated the value, both source and key are None. 971 972 Parameters 973 ---------- 974 constraint: Constraint, which must be in the subtree of a transformed 975 Disjunct 976 """ 977 transBlock = _get_constraint_transBlock(constraint) 978 # This is a KeyError if it fails, but it is also my fault if it 979 # fails... (That is, it's a bug in the mapping.) 980 return transBlock.bigm_src[constraint] 981 982 def get_M_value(self, constraint): 983 """Returns the M values used to transform constraint. Return is a tuple: 984 (lower_M_value, upper_M_value). Either can be None if constraint does 985 not have a lower or upper bound, respectively. 986 987 Parameters 988 ---------- 989 constraint: Constraint, which must be in the subtree of a transformed 990 Disjunct 991 """ 992 transBlock = _get_constraint_transBlock(constraint) 993 # This is a KeyError if it fails, but it is also my fault if it 994 # fails... (That is, it's a bug in the mapping.) 995 lower, upper = transBlock.bigm_src[constraint] 996 return (lower[0], upper[0]) 997