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 11import logging 12 13import pyomo.common.config as cfg 14from pyomo.common import deprecated 15from pyomo.common.collections import ComponentMap, ComponentSet 16from pyomo.common.log import is_debug_set 17from pyomo.common.modeling import unique_component_name 18from pyomo.core.expr.numvalue import ZeroConstant 19import pyomo.core.expr.current as EXPR 20from pyomo.core.base import Transformation, TransformationFactory, Reference 21from pyomo.core import ( 22 Block, BooleanVar, Connector, Constraint, Param, Set, SetOf, Suffix, Var, 23 Expression, SortComponents, TraversalStrategy, 24 Any, RangeSet, Reals, value, NonNegativeIntegers, LogicalConstraint, 25) 26from pyomo.gdp import Disjunct, Disjunction, GDP_Error 27from pyomo.gdp.util import ( _warn_for_active_logical_constraint, 28 clone_without_expression_components, target_list, 29 is_child_of, get_src_disjunction, 30 get_src_constraint, get_transformed_constraints, 31 get_src_disjunct, _warn_for_active_disjunction, 32 _warn_for_active_disjunct, preprocess_targets) 33from pyomo.network import Port 34from functools import wraps 35from weakref import ref as weakref_ref 36 37logger = logging.getLogger('pyomo.gdp.hull') 38 39NAME_BUFFER = {} 40 41@TransformationFactory.register( 42 'gdp.hull', 43 doc="Relax disjunctive model by forming the hull reformulation.") 44class Hull_Reformulation(Transformation): 45 """Relax disjunctive model by forming the hull reformulation. 46 47 Relaxes a disjunctive model into an algebraic model by forming the 48 hull reformulation of each disjunction. 49 50 This transformation accepts the following keyword arguments: 51 52 Parameters 53 ---------- 54 perspective_function : str 55 The perspective function used for the disaggregated variables. 56 Must be one of 'FurmanSawayaGrossmann' (default), 57 'LeeGrossmann', or 'GrossmannLee' 58 EPS : float 59 The value to use for epsilon [default: 1e-4] 60 targets : (block, disjunction, or list of those types) 61 The targets to transform. This can be a block, disjunction, or a 62 list of blocks and Disjunctions [default: the instance] 63 64 The transformation will create a new Block with a unique 65 name beginning "_pyomo_gdp_hull_reformulation". 66 The block will have a dictionary "_disaggregatedVarMap: 67 'srcVar': ComponentMap(<src var>:<disaggregated var>), 68 'disaggregatedVar': ComponentMap(<disaggregated var>:<src var>) 69 70 It will also have a ComponentMap "_bigMConstraintMap": 71 72 <disaggregated var>:<bounds constraint> 73 74 Last, it will contain an indexed Block named "relaxedDisjuncts", 75 which will hold the relaxed disjuncts. This block is indexed by 76 an integer indicating the order in which the disjuncts were relaxed. 77 Each block has a dictionary "_constraintMap": 78 79 'srcConstraints': ComponentMap(<transformed constraint>: 80 <src constraint>), 81 'transformedConstraints':ComponentMap(<src constraint container>: 82 <transformed constraint container>, 83 <src constraintData>: 84 [<transformed constraintDatas>] 85 ) 86 87 All transformed Disjuncts will have a pointer to the block their transformed 88 constraints are on, and all transformed Disjunctions will have a 89 pointer to the corresponding OR or XOR constraint. 90 91 The _pyomo_gdp_hull_reformulation block will have a ComponentMap 92 "_disaggregationConstraintMap": 93 <src var>:ComponentMap(<srcDisjunction>: <disaggregation constraint>) 94 95 """ 96 97 98 CONFIG = cfg.ConfigBlock('gdp.hull') 99 CONFIG.declare('targets', cfg.ConfigValue( 100 default=None, 101 domain=target_list, 102 description="target or list of targets that will be relaxed", 103 doc=""" 104 105 This specifies the target or list of targets to relax as either a 106 component or a list of components. If None (default), the entire model 107 is transformed. Note that if the transformation is done out of place, 108 the list of targets should be attached to the model before it is cloned, 109 and the list will specify the targets on the cloned instance.""" 110 )) 111 CONFIG.declare('perspective function', cfg.ConfigValue( 112 default='FurmanSawayaGrossmann', 113 domain=cfg.In(['FurmanSawayaGrossmann','LeeGrossmann','GrossmannLee']), 114 description='perspective function used for variable disaggregation', 115 doc=""" 116 The perspective function used for variable disaggregation 117 118 "LeeGrossmann" is the original NL convex hull from Lee & 119 Grossmann (2000) [1]_, which substitutes nonlinear constraints 120 121 h_ik(x) <= 0 122 123 with 124 125 x_k = sum( nu_ik ) 126 y_ik * h_ik( nu_ik/y_ik ) <= 0 127 128 "GrossmannLee" is an updated formulation from Grossmann & 129 Lee (2003) [2]_, which avoids divide-by-0 errors by using: 130 131 x_k = sum( nu_ik ) 132 (y_ik + eps) * h_ik( nu_ik/(y_ik + eps) ) <= 0 133 134 "FurmanSawayaGrossmann" (default) is an improved relaxation [3]_ 135 that is exact at 0 and 1 while avoiding numerical issues from 136 the Lee & Grossmann formulation by using: 137 138 x_k = sum( nu_ik ) 139 ((1-eps)*y_ik + eps) * h_ik( nu_ik/((1-eps)*y_ik + eps) ) \ 140 - eps * h_ki(0) * ( 1-y_ik ) <= 0 141 142 References 143 ---------- 144 .. [1] Lee, S., & Grossmann, I. E. (2000). New algorithms for 145 nonlinear generalized disjunctive programming. Computers and 146 Chemical Engineering, 24, 2125-2141 147 148 .. [2] Grossmann, I. E., & Lee, S. (2003). Generalized disjunctive 149 programming: Nonlinear convex hull relaxation and algorithms. 150 Computational Optimization and Applications, 26, 83-100. 151 152 .. [3] Furman, K., Sawaya, N., and Grossmann, I. A computationally 153 useful algebraic representation of nonlinear disjunctive convex 154 sets using the perspective function. Optimization Online 155 (2016). http://www.optimization-online.org/DB_HTML/2016/07/5544.html. 156 """ 157 )) 158 CONFIG.declare('EPS', cfg.ConfigValue( 159 default=1e-4, 160 domain=cfg.PositiveFloat, 161 description="Epsilon value to use in perspective function", 162 )) 163 CONFIG.declare('assume_fixed_vars_permanent', cfg.ConfigValue( 164 default=False, 165 domain=bool, 166 description="Boolean indicating whether or not to transform so that the " 167 "the transformed model will still be valid when fixed Vars are unfixed.", 168 doc=""" 169 If True, the transformation will not disaggregate fixed variables. 170 This means that if a fixed variable is unfixed after transformation, 171 the transformed model is no longer valid. By default, the transformation 172 will disagregate fixed variables so that any later fixing and unfixing 173 will be valid in the transformed model. 174 """ 175 )) 176 177 def __init__(self): 178 super(Hull_Reformulation, self).__init__() 179 self.handlers = { 180 Constraint : self._transform_constraint, 181 Var : False, 182 BooleanVar: False, 183 Connector : False, 184 Expression : False, 185 Param : False, 186 Set : False, 187 SetOf : False, 188 RangeSet: False, 189 Suffix : False, 190 Disjunction: self._warn_for_active_disjunction, 191 Disjunct: self._warn_for_active_disjunct, 192 Block: self._transform_block_on_disjunct, 193 LogicalConstraint: self._warn_for_active_logical_statement, 194 Port: False, 195 } 196 self._generate_debug_messages = False 197 198 def _add_local_vars(self, block, local_var_dict): 199 localVars = block.component('LocalVars') 200 if type(localVars) is Suffix: 201 for disj, var_list in localVars.items(): 202 if local_var_dict.get(disj) is None: 203 local_var_dict[disj] = ComponentSet(var_list) 204 else: 205 local_var_dict[disj].update(var_list) 206 207 def _get_local_var_suffixes(self, block, local_var_dict): 208 # You can specify suffixes on any block (disjuncts included). This method 209 # starts from a Disjunct (presumably) and checks for a LocalVar suffixes 210 # going both up and down the tree, adding them into the dictionary that 211 # is the second argument. 212 213 # first look beneath where we are (there could be Blocks on this 214 # disjunct) 215 for b in block.component_data_objects(Block, descend_into=(Block), 216 active=True, 217 sort=SortComponents.deterministic): 218 self._add_local_vars(b, local_var_dict) 219 # now traverse upwards and get what's above 220 while block is not None: 221 self._add_local_vars(block, local_var_dict) 222 block = block.parent_block() 223 224 return local_var_dict 225 226 def _apply_to(self, instance, **kwds): 227 assert not NAME_BUFFER 228 try: 229 self._apply_to_impl(instance, **kwds) 230 finally: 231 # Clear the global name buffer now that we are done 232 NAME_BUFFER.clear() 233 234 def _apply_to_impl(self, instance, **kwds): 235 if not instance.ctype in (Block, Disjunct): 236 raise GDP_Error("Transformation called on %s of type %s. 'instance' " 237 "must be a ConcreteModel, Block, or Disjunct (in " 238 "the case of nested disjunctions)." % 239 (instance.name, instance.ctype)) 240 241 self._config = self.CONFIG(kwds.pop('options', {})) 242 self._config.set_value(kwds) 243 self._generate_debug_messages = is_debug_set(logger) 244 245 targets = self._config.targets 246 if targets is None: 247 targets = ( instance, ) 248 else: 249 # we need to preprocess targets to make sure that if there are any 250 # disjunctions in targets that their disjuncts appear before them in 251 # the list. 252 targets = preprocess_targets(targets) 253 knownBlocks = {} 254 for t in targets: 255 # check that t is in fact a child of instance 256 if not is_child_of(parent=instance, child=t, 257 knownBlocks=knownBlocks): 258 raise GDP_Error( 259 "Target '%s' is not a component on instance '%s'!" 260 % (t.name, instance.name)) 261 elif t.ctype is Disjunction: 262 if t.is_indexed(): 263 self._transform_disjunction(t) 264 else: 265 self._transform_disjunctionData(t, t.index()) 266 elif t.ctype in (Block, Disjunct): 267 if t.is_indexed(): 268 self._transform_block(t) 269 else: 270 self._transform_blockData(t) 271 else: 272 raise GDP_Error( 273 "Target '%s' was not a Block, Disjunct, or Disjunction. " 274 "It was of type %s and can't be transformed." 275 % (t.name, type(t)) ) 276 277 def _add_transformation_block(self, instance): 278 # make a transformation block on instance where we will store 279 # transformed components 280 transBlockName = unique_component_name( 281 instance, 282 '_pyomo_gdp_hull_reformulation') 283 transBlock = Block() 284 instance.add_component(transBlockName, transBlock) 285 transBlock.relaxedDisjuncts = Block(NonNegativeIntegers) 286 transBlock.lbub = Set(initialize = ['lb','ub','eq']) 287 # Map between disaggregated variables and their 288 # originals 289 transBlock._disaggregatedVarMap = { 290 'srcVar': ComponentMap(), 291 'disaggregatedVar': ComponentMap(), 292 } 293 # Map between disaggregated variables and their lb*indicator <= var <= 294 # ub*indicator constraints 295 transBlock._bigMConstraintMap = ComponentMap() 296 # We will store all of the disaggregation constraints for any 297 # Disjunctions we transform onto this block here. 298 transBlock.disaggregationConstraints = Constraint(NonNegativeIntegers, 299 Any) 300 301 # This will map from srcVar to a map of srcDisjunction to the 302 # disaggregation constraint corresponding to srcDisjunction 303 transBlock._disaggregationConstraintMap = ComponentMap() 304 305 # we are going to store some of the disaggregated vars directly here 306 # when we have vars that don't appear in every disjunct 307 transBlock._disaggregatedVars = Var(NonNegativeIntegers, dense=False) 308 transBlock._boundsConstraints = Constraint(NonNegativeIntegers, 309 transBlock.lbub) 310 311 return transBlock 312 313 def _transform_block(self, obj): 314 for i in sorted(obj.keys()): 315 self._transform_blockData(obj[i]) 316 317 def _transform_blockData(self, obj): 318 # Transform every (active) disjunction in the block 319 for disjunction in obj.component_objects( 320 Disjunction, 321 active=True, 322 sort=SortComponents.deterministic, 323 descend_into=(Block,Disjunct), 324 descent_order=TraversalStrategy.PostfixDFS): 325 self._transform_disjunction(disjunction) 326 327 def _add_xor_constraint(self, disjunction, transBlock): 328 # Put XOR constraint on the transformation block 329 330 # We never do this for just a DisjunctionData because we need 331 # to know about the index set of its parent component. So if 332 # we called this on a DisjunctionData, we did something wrong. 333 assert isinstance(disjunction, Disjunction) 334 335 # check if the constraint already exists 336 if disjunction._algebraic_constraint is not None: 337 return disjunction._algebraic_constraint() 338 339 # add the XOR (or OR) constraints to parent block (with 340 # unique name) It's indexed if this is an 341 # IndexedDisjunction, not otherwise 342 orC = Constraint(disjunction.index_set()) 343 transBlock.add_component( 344 unique_component_name(transBlock, 345 disjunction.getname(fully_qualified=True, 346 name_buffer=NAME_BUFFER) +\ 347 '_xor'), orC) 348 disjunction._algebraic_constraint = weakref_ref(orC) 349 350 return orC 351 352 def _transform_disjunction(self, obj): 353 # NOTE: this check is actually necessary because it's possible we go 354 # straight to this function when we use targets. 355 if not obj.active: 356 return 357 358 # put the transformation block on the parent block of the Disjunction, 359 # unless this is a disjunction we have seen in a prior call to hull, in 360 # which case we will use the same transformation block we created 361 # before. 362 if obj._algebraic_constraint is not None: 363 transBlock = obj._algebraic_constraint().parent_block() 364 else: 365 transBlock = self._add_transformation_block(obj.parent_block()) 366 # and create the xor constraint 367 xorConstraint = self._add_xor_constraint(obj, transBlock) 368 369 # create the disjunction constraint and disaggregation 370 # constraints and then relax each of the disjunctionDatas 371 for i in sorted(obj.keys()): 372 self._transform_disjunctionData(obj[i], i, transBlock) 373 374 # deactivate so the writers will be happy 375 obj.deactivate() 376 377 def _transform_disjunctionData(self, obj, index, transBlock=None): 378 if not obj.active: 379 return 380 # Hull reformulation doesn't work if this is an OR constraint. So if 381 # xor is false, give up 382 if not obj.xor: 383 raise GDP_Error("Cannot do hull reformulation for " 384 "Disjunction '%s' with OR constraint. " 385 "Must be an XOR!" % obj.name) 386 387 if transBlock is None: 388 # It's possible that we have already created a transformation block 389 # for another disjunctionData from this same container. If that's 390 # the case, let's use the same transformation block. (Else it will 391 # be really confusing that the XOR constraint goes to that old block 392 # but we create a new one here.) 393 if obj.parent_component()._algebraic_constraint is not None: 394 transBlock = obj.parent_component()._algebraic_constraint().\ 395 parent_block() 396 else: 397 transBlock = self._add_transformation_block(obj.parent_block()) 398 399 parent_component = obj.parent_component() 400 401 orConstraint = self._add_xor_constraint(parent_component, transBlock) 402 disaggregationConstraint = transBlock.disaggregationConstraints 403 disaggregationConstraintMap = transBlock._disaggregationConstraintMap 404 disaggregatedVars = transBlock._disaggregatedVars 405 disaggregated_var_bounds = transBlock._boundsConstraints 406 407 # Just because it's unlikely this is what someone meant to do... 408 if len(obj.disjuncts) == 0: 409 raise GDP_Error("Disjunction '%s' is empty. This is " 410 "likely indicative of a modeling error." % 411 obj.getname(fully_qualified=True, 412 name_buffer=NAME_BUFFER)) 413 414 # We first go through and collect all the variables that we 415 # are going to disaggregate. 416 varOrder_set = ComponentSet() 417 varOrder = [] 418 varsByDisjunct = ComponentMap() 419 localVarsByDisjunct = ComponentMap() 420 include_fixed_vars = not self._config.assume_fixed_vars_permanent 421 for disjunct in obj.disjuncts: 422 disjunctVars = varsByDisjunct[disjunct] = ComponentSet() 423 # create the key for each disjunct now 424 transBlock._disaggregatedVarMap['disaggregatedVar'][ 425 disjunct] = ComponentMap() 426 for cons in disjunct.component_data_objects( 427 Constraint, 428 active = True, 429 sort=SortComponents.deterministic, 430 descend_into=Block): 431 # [ESJ 02/14/2020] By default, we disaggregate fixed variables 432 # on the philosophy that fixing is not a promise for the future 433 # and we are mathematically wrong if we don't transform these 434 # correctly and someone later unfixes them and keeps playing 435 # with their transformed model. However, the user may have set 436 # assume_fixed_vars_permanent to True in which case we will skip 437 # them 438 for var in EXPR.identify_variables( 439 cons.body, include_fixed=include_fixed_vars): 440 # Note the use of a list so that we will 441 # eventually disaggregate the vars in a 442 # deterministic order (the order that we found 443 # them) 444 disjunctVars.add(var) 445 if not var in varOrder_set: 446 varOrder.append(var) 447 varOrder_set.add(var) 448 449 # check for LocalVars Suffix 450 localVarsByDisjunct = self._get_local_var_suffixes( 451 disjunct, localVarsByDisjunct) 452 453 # We will disaggregate all variables which are not explicitly declared 454 # as being local. Note however, that we do declare our own disaggregated 455 # variables as local, so they will not be re-disaggregated. 456 varSet = [] 457 varSet = {disj: [] for disj in obj.disjuncts} 458 # Note that variables are local with respect to a Disjunct. We deal with 459 # them here to do some error checking (if something is obviously not 460 # local since it is used in multiple Disjuncts in this Disjunction) and 461 # also to get a deterministic order in which to process them when we 462 # transform the Disjuncts: Values of localVarsByDisjunct are 463 # ComponentSets, so we need this for determinism (we iterate through the 464 # localVars of a Disjunct later) 465 localVars = ComponentMap() 466 varsToDisaggregate = [] 467 disjunctsVarAppearsIn = ComponentMap() 468 for var in varOrder: 469 disjuncts = disjunctsVarAppearsIn[var] = [d for d in varsByDisjunct 470 if var in 471 varsByDisjunct[d]] 472 # clearly not local if used in more than one disjunct 473 if len(disjuncts) > 1: 474 if self._generate_debug_messages: 475 logger.debug("Assuming '%s' is not a local var since it is" 476 "used in multiple disjuncts." % 477 var.getname(fully_qualified=True, 478 name_buffer=NAME_BUFFER)) 479 for disj in disjuncts: 480 varSet[disj].append(var) 481 varsToDisaggregate.append(var) 482 # disjuncts is a list of length 1 483 elif localVarsByDisjunct.get(disjuncts[0]) is not None: 484 if var in localVarsByDisjunct[disjuncts[0]]: 485 localVars_thisDisjunct = localVars.get(disjuncts[0]) 486 if localVars_thisDisjunct is not None: 487 localVars[disjuncts[0]].append(var) 488 else: 489 localVars[disjuncts[0]] = [var] 490 else: 491 # It's not local to this Disjunct 492 varSet[disjuncts[0]].append(var) 493 varsToDisaggregate.append(var) 494 else: 495 # We don't even have have any local vars for this Disjunct. 496 varSet[disjuncts[0]].append(var) 497 varsToDisaggregate.append(var) 498 499 # Now that we know who we need to disaggregate, we will do it 500 # while we also transform the disjuncts. 501 local_var_set = self._get_local_var_set(obj) 502 or_expr = 0 503 for disjunct in obj.disjuncts: 504 or_expr += disjunct.indicator_var.get_associated_binary() 505 self._transform_disjunct(disjunct, transBlock, varSet[disjunct], 506 localVars.get(disjunct, []), local_var_set) 507 orConstraint.add(index, (or_expr, 1)) 508 # map the DisjunctionData to its XOR constraint to mark it as 509 # transformed 510 obj._algebraic_constraint = weakref_ref(orConstraint[index]) 511 512 # add the reaggregation constraints 513 for i, var in enumerate(varsToDisaggregate): 514 # There are two cases here: Either the var appeared in every 515 # disjunct in the disjunction, or it didn't. If it did, there's 516 # nothing special to do: All of the disaggregated variables have 517 # been created, and we can just proceed and make this constraint. If 518 # it didn't, we need one more disaggregated variable, correctly 519 # defined. And then we can make the constraint. 520 if len(disjunctsVarAppearsIn[var]) < len(obj.disjuncts): 521 # create one more disaggregated var 522 idx = len(disaggregatedVars) 523 disaggregated_var = disaggregatedVars[idx] 524 # mark this as local because we won't re-disaggregate if this is 525 # a nested disjunction 526 if local_var_set is not None: 527 local_var_set.append(disaggregatedVar) 528 var_free = 1 - sum(disj.indicator_var.get_associated_binary() 529 for disj in disjunctsVarAppearsIn[var]) 530 self._declare_disaggregated_var_bounds(var, disaggregated_var, 531 obj, 532 disaggregated_var_bounds, 533 (idx,'lb'), (idx,'ub'), 534 var_free) 535 # maintain the mappings 536 for disj in obj.disjuncts: 537 # Because we called _transform_disjunct above, we know that 538 # if this isn't transformed it is because it was cleanly 539 # deactivated, and we can just skip it. 540 if disj._transformation_block is not None and \ 541 disj not in disjunctsVarAppearsIn[var]: 542 relaxationBlock = disj._transformation_block().\ 543 parent_block() 544 relaxationBlock._bigMConstraintMap[ 545 disaggregated_var] = Reference( 546 disaggregated_var_bounds[idx, :]) 547 relaxationBlock._disaggregatedVarMap['srcVar'][ 548 disaggregated_var] = var 549 relaxationBlock._disaggregatedVarMap[ 550 'disaggregatedVar'][disj][var] = disaggregated_var 551 552 disaggregatedExpr = disaggregated_var 553 else: 554 disaggregatedExpr = 0 555 for disjunct in disjunctsVarAppearsIn[var]: 556 if disjunct._transformation_block is None: 557 # Because we called _transform_disjunct above, we know that 558 # if this isn't transformed it is because it was cleanly 559 # deactivated, and we can just skip it. 560 continue 561 562 disaggregatedVar = disjunct._transformation_block().\ 563 parent_block()._disaggregatedVarMap[ 564 'disaggregatedVar'][disjunct][var] 565 disaggregatedExpr += disaggregatedVar 566 567 disaggregationConstraint.add((i, index), var == disaggregatedExpr) 568 # and update the map so that we can find this later. We index by 569 # variable and the particular disjunction because there is a 570 # different one for each disjunction 571 if disaggregationConstraintMap.get(var) is not None: 572 disaggregationConstraintMap[var][obj] = disaggregationConstraint[ 573 (i, index)] 574 else: 575 thismap = disaggregationConstraintMap[var] = ComponentMap() 576 thismap[obj] = disaggregationConstraint[(i, index)] 577 578 # deactivate for the writers 579 obj.deactivate() 580 581 def _transform_disjunct(self, obj, transBlock, varSet, localVars, 582 local_var_set): 583 # deactivated should only come from the user 584 if not obj.active: 585 if obj.indicator_var.is_fixed(): 586 if not value(obj.indicator_var): 587 # The user cleanly deactivated the disjunct: there 588 # is nothing for us to do here. 589 return 590 else: 591 raise GDP_Error( 592 "The disjunct '%s' is deactivated, but the " 593 "indicator_var is fixed to %s. This makes no sense." 594 % ( obj.name, value(obj.indicator_var) )) 595 if obj._transformation_block is None: 596 raise GDP_Error( 597 "The disjunct '%s' is deactivated, but the " 598 "indicator_var is not fixed and the disjunct does not " 599 "appear to have been relaxed. This makes no sense. " 600 "(If the intent is to deactivate the disjunct, fix its " 601 "indicator_var to False.)" 602 % ( obj.name, )) 603 604 if obj._transformation_block is not None: 605 # we've transformed it, which means this is the second time it's 606 # appearing in a Disjunction 607 raise GDP_Error( 608 "The disjunct '%s' has been transformed, but a disjunction " 609 "it appears in has not. Putting the same disjunct in " 610 "multiple disjunctions is not supported." % obj.name) 611 612 # create a relaxation block for this disjunct 613 relaxedDisjuncts = transBlock.relaxedDisjuncts 614 relaxationBlock = relaxedDisjuncts[len(relaxedDisjuncts)] 615 transBlock = relaxationBlock.parent_block() 616 617 relaxationBlock.localVarReferences = Block() 618 619 # Put the disaggregated variables all on their own block so that we can 620 # isolate the name collisions and still have complete control over the 621 # names on this block. (This is for peace of mind now, but will matter 622 # in the future for adding the binaries corresponding to Boolean 623 # indicator vars.) 624 relaxationBlock.disaggregatedVars = Block() 625 626 # add the map that will link back and forth between transformed 627 # constraints and their originals. 628 relaxationBlock._constraintMap = { 629 'srcConstraints': ComponentMap(), 630 'transformedConstraints': ComponentMap() 631 } 632 633 # add mappings to source disjunct (so we'll know we've relaxed) 634 obj._transformation_block = weakref_ref(relaxationBlock) 635 relaxationBlock._srcDisjunct = weakref_ref(obj) 636 637 # add the disaggregated variables and their bigm constraints 638 # to the relaxationBlock 639 for var in varSet: 640 disaggregatedVar = Var(within=Reals, initialize=var.value) 641 # naming conflicts are possible here since this is a bunch 642 # of variables from different blocks coming together, so we 643 # get a unique name 644 disaggregatedVarName = unique_component_name( 645 relaxationBlock.disaggregatedVars, 646 var.getname(fully_qualified=False, name_buffer=NAME_BUFFER)) 647 relaxationBlock.disaggregatedVars.add_component( 648 disaggregatedVarName, disaggregatedVar) 649 # mark this as local because we won't re-disaggregate if this is a 650 # nested disjunction 651 if local_var_set is not None: 652 local_var_set.append(disaggregatedVar) 653 654 # add the bigm constraint 655 bigmConstraint = Constraint(transBlock.lbub) 656 relaxationBlock.add_component( 657 disaggregatedVarName + "_bounds", bigmConstraint) 658 659 self._declare_disaggregated_var_bounds( 660 var, disaggregatedVar, obj, 661 bigmConstraint, 'lb', 'ub', 662 obj.indicator_var.get_associated_binary(), transBlock) 663 664 for var in localVars: 665 # we don't need to disaggregated, we can use this Var, but we do 666 # need to set up its bounds constraints. 667 668 # naming conflicts are possible here since this is a bunch 669 # of variables from different blocks coming together, so we 670 # get a unique name 671 conName = unique_component_name( 672 relaxationBlock, 673 var.getname(fully_qualified=False, name_buffer=NAME_BUFFER) + \ 674 "_bounds") 675 bigmConstraint = Constraint(transBlock.lbub) 676 relaxationBlock.add_component(conName, bigmConstraint) 677 678 self._declare_disaggregated_var_bounds( 679 var, var, obj, 680 bigmConstraint, 'lb', 'ub', 681 obj.indicator_var.get_associated_binary(), transBlock) 682 683 var_substitute_map = dict((id(v), newV) for v, newV in 684 transBlock._disaggregatedVarMap[ 685 'disaggregatedVar'][obj].items()) 686 zero_substitute_map = dict((id(v), ZeroConstant) for v, newV in \ 687 transBlock._disaggregatedVarMap[ 688 'disaggregatedVar'][obj].items()) 689 zero_substitute_map.update((id(v), ZeroConstant) for v in localVars) 690 691 # Transform each component within this disjunct 692 self._transform_block_components(obj, obj, var_substitute_map, 693 zero_substitute_map) 694 695 # deactivate disjunct so writers can be happy 696 obj._deactivate_without_fixing_indicator() 697 698 def _transform_block_components( self, block, disjunct, var_substitute_map, 699 zero_substitute_map): 700 # As opposed to bigm, in hull the only special thing we need to do for 701 # nested Disjunctions is to make sure that we move up local var 702 # references and also references to the disaggregated variables so that 703 # all will be accessible after we transform this Disjunct. The indicator 704 # variables and disaggregated variables of the inner disjunction will 705 # need to be disaggregated again, but the transformed constraints will 706 # not be. But this way nothing will get double-bigm-ed. (If an 707 # untransformed disjunction is lurking here, we will catch it below). 708 709 disjunctBlock = disjunct._transformation_block() 710 destinationBlock = disjunctBlock.parent_block() 711 for obj in block.component_data_objects( 712 Disjunction, 713 sort=SortComponents.deterministic, 714 descend_into=(Block)): 715 if obj.algebraic_constraint is None: 716 # This could be bad if it's active since that means its 717 # untransformed, but we'll wait to yell until the next loop 718 continue 719 # get this disjunction's relaxation block. 720 transBlock = obj.algebraic_constraint().parent_block() 721 722 self._transfer_var_references(transBlock, destinationBlock) 723 724 # add references to all local variables on block (including the 725 # indicator_var). Note that we do this after we have moved up the 726 # transformation blocks for nested disjunctions, so that we don't have 727 # duplicate references. 728 varRefBlock = disjunctBlock.localVarReferences 729 for v in block.component_objects(Var, descend_into=Block, active=None): 730 if len(v) > 0: 731 varRefBlock.add_component(unique_component_name( 732 varRefBlock, v.getname(fully_qualified=True, 733 name_buffer=NAME_BUFFER)), 734 Reference(v)) 735 736 # Look through the component map of block and transform everything we 737 # have a handler for. Yell if we don't know how to handle it. (Note that 738 # because we only iterate through active components, this means 739 # non-ActiveComponent types cannot have handlers.) 740 for obj in block.component_objects(active=True, descend_into=False): 741 handler = self.handlers.get(obj.ctype, None) 742 if not handler: 743 if handler is None: 744 raise GDP_Error( 745 "No hull transformation handler registered " 746 "for modeling components of type %s. If your " 747 "disjuncts contain non-GDP Pyomo components that " 748 "require transformation, please transform them first." 749 % obj.ctype ) 750 continue 751 # obj is what we are transforming, we pass disjunct 752 # through so that we will have access to the indicator 753 # variables down the line. 754 handler(obj, disjunct, var_substitute_map, zero_substitute_map) 755 756 def _declare_disaggregated_var_bounds(self, original_var, disaggregatedVar, 757 disjunct, bigmConstraint, lb_idx, 758 ub_idx, var_free_indicator, 759 transBlock=None): 760 # If transBlock is None then this is a disaggregated variable for 761 # multiple Disjuncts and we will handle the mappings separately. 762 lb = original_var.lb 763 ub = original_var.ub 764 if lb is None or ub is None: 765 raise GDP_Error("Variables that appear in disjuncts must be " 766 "bounded in order to use the hull " 767 "transformation! Missing bound for %s." 768 % (original_var.name)) 769 770 disaggregatedVar.setlb(min(0, lb)) 771 disaggregatedVar.setub(max(0, ub)) 772 773 if lb: 774 bigmConstraint.add( 775 lb_idx, var_free_indicator*lb <= disaggregatedVar) 776 if ub: 777 bigmConstraint.add( 778 ub_idx, disaggregatedVar <= ub*var_free_indicator) 779 780 # store the mappings from variables to their disaggregated selves on 781 # the transformation block. 782 if transBlock is not None: 783 transBlock._disaggregatedVarMap['disaggregatedVar'][disjunct][ 784 original_var] = disaggregatedVar 785 transBlock._disaggregatedVarMap['srcVar'][ 786 disaggregatedVar] = original_var 787 transBlock._bigMConstraintMap[disaggregatedVar] = bigmConstraint 788 789 def _get_local_var_set(self, disjunction): 790 # add Suffix to the relaxation block that disaggregated variables are 791 # local (in case this is nested in another Disjunct) 792 local_var_set = None 793 parent_disjunct = disjunction.parent_block() 794 while parent_disjunct is not None: 795 if parent_disjunct.ctype is Disjunct: 796 break 797 parent_disjunct = parent_disjunct.parent_block() 798 if parent_disjunct is not None: 799 # This limits the cases that a user is allowed to name something 800 # (other than a Suffix) 'LocalVars' on a Disjunct. But I am assuming 801 # that the Suffix has to be somewhere above the disjunct in the 802 # tree, so I can't put it on a Block that I own. And if I'm coopting 803 # something of theirs, it may as well be here. 804 self._add_local_var_suffix(parent_disjunct) 805 if parent_disjunct.LocalVars.get(parent_disjunct) is None: 806 parent_disjunct.LocalVars[parent_disjunct] = [] 807 local_var_set = parent_disjunct.LocalVars[parent_disjunct] 808 809 return local_var_set 810 811 def _transfer_var_references(self, fromBlock, toBlock): 812 disjunctList = toBlock.relaxedDisjuncts 813 for idx, disjunctBlock in fromBlock.relaxedDisjuncts.items(): 814 # move all the of the local var references 815 newblock = disjunctList[len(disjunctList)] 816 newblock.localVarReferences = Block() 817 newblock.localVarReferences.transfer_attributes_from( 818 disjunctBlock.localVarReferences) 819 820 def _warn_for_active_disjunction( self, disjunction, disjunct, 821 var_substitute_map, zero_substitute_map): 822 _warn_for_active_disjunction(disjunction, disjunct, NAME_BUFFER) 823 824 def _warn_for_active_disjunct( self, innerdisjunct, outerdisjunct, 825 var_substitute_map, zero_substitute_map): 826 _warn_for_active_disjunct(innerdisjunct, outerdisjunct, NAME_BUFFER) 827 828 def _warn_for_active_logical_statement( 829 self, logical_statment, disjunct, var_substitute_map, 830 zero_substitute_map): 831 _warn_for_active_logical_constraint(logical_statment, disjunct, 832 NAME_BUFFER) 833 834 def _transform_block_on_disjunct( self, block, disjunct, var_substitute_map, 835 zero_substitute_map): 836 # We look through everything on the component map of the block 837 # and transform it just as we would if it was on the disjunct 838 # directly. (We are passing the disjunct through so that when 839 # we find constraints, _transform_constraint will have access to 840 # the correct indicator variable. 841 for i in sorted(block.keys()): 842 self._transform_block_components( block[i], disjunct, 843 var_substitute_map, 844 zero_substitute_map) 845 846 def _transform_constraint(self, obj, disjunct, var_substitute_map, 847 zero_substitute_map): 848 # we will put a new transformed constraint on the relaxation block. 849 relaxationBlock = disjunct._transformation_block() 850 transBlock = relaxationBlock.parent_block() 851 constraintMap = relaxationBlock._constraintMap 852 853 # Though rare, it is possible to get naming conflicts here 854 # since constraints from all blocks are getting moved onto the 855 # same block. So we get a unique name 856 name = unique_component_name(relaxationBlock, obj.getname( 857 fully_qualified=True, name_buffer=NAME_BUFFER)) 858 859 if obj.is_indexed(): 860 newConstraint = Constraint(obj.index_set(), transBlock.lbub) 861 else: 862 newConstraint = Constraint(transBlock.lbub) 863 relaxationBlock.add_component(name, newConstraint) 864 # map the containers: 865 # add mapping of original constraint to transformed constraint 866 if obj.is_indexed(): 867 constraintMap['transformedConstraints'][obj] = newConstraint 868 # add mapping of transformed constraint container back to original 869 # constraint container (or ScalarConstraint) 870 constraintMap['srcConstraints'][newConstraint] = obj 871 872 for i in sorted(obj.keys()): 873 c = obj[i] 874 if not c.active: 875 continue 876 877 NL = c.body.polynomial_degree() not in (0,1) 878 EPS = self._config.EPS 879 mode = self._config.perspective_function 880 881 # We need to evaluate the expression at the origin *before* 882 # we substitute the expression variables with the 883 # disaggregated variables 884 if not NL or mode == "FurmanSawayaGrossmann": 885 h_0 = clone_without_expression_components( 886 c.body, substitute=zero_substitute_map) 887 888 y = disjunct.binary_indicator_var 889 if NL: 890 if mode == "LeeGrossmann": 891 sub_expr = clone_without_expression_components( 892 c.body, 893 substitute=dict( 894 (var, subs/y) 895 for var, subs in var_substitute_map.items() ) 896 ) 897 expr = sub_expr * y 898 elif mode == "GrossmannLee": 899 sub_expr = clone_without_expression_components( 900 c.body, 901 substitute=dict( 902 (var, subs/(y + EPS)) 903 for var, subs in var_substitute_map.items() ) 904 ) 905 expr = (y + EPS) * sub_expr 906 elif mode == "FurmanSawayaGrossmann": 907 sub_expr = clone_without_expression_components( 908 c.body, 909 substitute=dict( 910 (var, subs/((1 - EPS)*y + EPS)) 911 for var, subs in var_substitute_map.items() ) 912 ) 913 expr = ((1-EPS)*y + EPS)*sub_expr - EPS*h_0*(1-y) 914 else: 915 raise RuntimeError("Unknown NL Hull mode") 916 else: 917 expr = clone_without_expression_components( 918 c.body, substitute=var_substitute_map) 919 920 if c.equality: 921 if NL: 922 # ESJ TODO: This can't happen right? This is the only 923 # obvious case where someone has messed up, but this has to 924 # be nonconvex, right? Shouldn't we tell them? 925 newConsExpr = expr == c.lower*y 926 else: 927 v = list(EXPR.identify_variables(expr)) 928 if len(v) == 1 and not c.lower: 929 # Setting a variable to 0 in a disjunct is 930 # *very* common. We should recognize that in 931 # that structure, the disaggregated variable 932 # will also be fixed to 0. 933 v[0].fix(0) 934 # ESJ: If you ask where the transformed constraint is, 935 # the answer is nowhere. Really, it is in the bounds of 936 # this variable, so I'm going to return 937 # it. Alternatively we could return an empty list, but I 938 # think I like this better. 939 constraintMap['transformedConstraints'][c] = [v[0]] 940 # Reverse map also (this is strange) 941 constraintMap['srcConstraints'][v[0]] = c 942 continue 943 newConsExpr = expr - (1-y)*h_0 == c.lower*y 944 945 if obj.is_indexed(): 946 newConstraint.add((i, 'eq'), newConsExpr) 947 # map the _ConstraintDatas (we mapped the container above) 948 constraintMap[ 949 'transformedConstraints'][c] = [newConstraint[i,'eq']] 950 constraintMap['srcConstraints'][newConstraint[i,'eq']] = c 951 else: 952 newConstraint.add('eq', newConsExpr) 953 # map to the _ConstraintData (And yes, for 954 # ScalarConstraints, this is overwriting the map to the 955 # container we made above, and that is what I want to 956 # happen. ScalarConstraints will map to lists. For 957 # IndexedConstraints, we can map the container to the 958 # container, but more importantly, we are mapping the 959 # _ConstraintDatas to each other above) 960 constraintMap[ 961 'transformedConstraints'][c] = [newConstraint['eq']] 962 constraintMap['srcConstraints'][newConstraint['eq']] = c 963 964 continue 965 966 if c.lower is not None: 967 if self._generate_debug_messages: 968 _name = c.getname( 969 fully_qualified=True, name_buffer=NAME_BUFFER) 970 logger.debug("GDP(Hull): Transforming constraint " + 971 "'%s'", _name) 972 if NL: 973 newConsExpr = expr >= c.lower*y 974 else: 975 newConsExpr = expr - (1-y)*h_0 >= c.lower*y 976 977 if obj.is_indexed(): 978 newConstraint.add((i, 'lb'), newConsExpr) 979 constraintMap[ 980 'transformedConstraints'][c] = [newConstraint[i,'lb']] 981 constraintMap['srcConstraints'][newConstraint[i,'lb']] = c 982 else: 983 newConstraint.add('lb', newConsExpr) 984 constraintMap[ 985 'transformedConstraints'][c] = [newConstraint['lb']] 986 constraintMap['srcConstraints'][newConstraint['lb']] = c 987 988 if c.upper is not None: 989 if self._generate_debug_messages: 990 _name = c.getname( 991 fully_qualified=True, name_buffer=NAME_BUFFER) 992 logger.debug("GDP(Hull): Transforming constraint " + 993 "'%s'", _name) 994 if NL: 995 newConsExpr = expr <= c.upper*y 996 else: 997 newConsExpr = expr - (1-y)*h_0 <= c.upper*y 998 999 if obj.is_indexed(): 1000 newConstraint.add((i, 'ub'), newConsExpr) 1001 # map (have to account for fact we might have created list 1002 # above 1003 transformed = constraintMap['transformedConstraints'].get(c) 1004 if transformed is not None: 1005 transformed.append(newConstraint[i,'ub']) 1006 else: 1007 constraintMap['transformedConstraints'][ 1008 c] = [newConstraint[i,'ub']] 1009 constraintMap['srcConstraints'][newConstraint[i,'ub']] = c 1010 else: 1011 newConstraint.add('ub', newConsExpr) 1012 transformed = constraintMap['transformedConstraints'].get(c) 1013 if transformed is not None: 1014 transformed.append(newConstraint['ub']) 1015 else: 1016 constraintMap['transformedConstraints'][ 1017 c] = [newConstraint['ub']] 1018 constraintMap['srcConstraints'][newConstraint['ub']] = c 1019 1020 # deactivate now that we have transformed 1021 obj.deactivate() 1022 1023 def _add_local_var_suffix(self, disjunct): 1024 # If the Suffix is there, we will borrow it. If not, we make it. If it's 1025 # something else, we complain. 1026 localSuffix = disjunct.component("LocalVars") 1027 if localSuffix is None: 1028 disjunct.LocalVars = Suffix(direction=Suffix.LOCAL) 1029 else: 1030 if localSuffix.ctype is Suffix: 1031 return 1032 raise GDP_Error("A component called 'LocalVars' is declared on " 1033 "Disjunct %s, but it is of type %s, not Suffix." 1034 % (disjunct.getname(fully_qualified=True, 1035 name_buffer=NAME_BUFFER), 1036 localSuffix.ctype)) 1037 1038 # These are all functions to retrieve transformed components from 1039 # original ones and vice versa. 1040 1041 @wraps(get_src_disjunct) 1042 def get_src_disjunct(self, transBlock): 1043 return get_src_disjunct(transBlock) 1044 1045 @wraps(get_src_disjunction) 1046 def get_src_disjunction(self, xor_constraint): 1047 return get_src_disjunction(xor_constraint) 1048 1049 @wraps(get_src_constraint) 1050 def get_src_constraint(self, transformedConstraint): 1051 return get_src_constraint(transformedConstraint) 1052 1053 @wraps(get_transformed_constraints) 1054 def get_transformed_constraints(self, srcConstraint): 1055 return get_transformed_constraints(srcConstraint) 1056 1057 def get_disaggregated_var(self, v, disjunct): 1058 """ 1059 Returns the disaggregated variable corresponding to the Var v and the 1060 Disjunct disjunct. 1061 1062 If v is a local variable, this method will return v. 1063 1064 Parameters 1065 ---------- 1066 v: a Var which appears in a constraint in a transformed Disjunct 1067 disjunct: a transformed Disjunct in which v appears 1068 """ 1069 if disjunct._transformation_block is None: 1070 raise GDP_Error("Disjunct '%s' has not been transformed" 1071 % disjunct.name) 1072 transBlock = disjunct._transformation_block().parent_block() 1073 try: 1074 return transBlock._disaggregatedVarMap['disaggregatedVar'][ 1075 disjunct][v] 1076 except: 1077 logger.error("It does not appear '%s' is a " 1078 "variable which appears in disjunct '%s'" 1079 % (v.name, disjunct.name)) 1080 raise 1081 1082 def get_src_var(self, disaggregated_var): 1083 """ 1084 Returns the original model variable to which disaggregated_var 1085 corresponds. 1086 1087 Parameters 1088 ---------- 1089 disaggregated_var: a Var which was created by the hull 1090 transformation as a disaggregated variable 1091 (and so appears on a transformation block 1092 of some Disjunct) 1093 """ 1094 msg = ("'%s' does not appear to be a " 1095 "disaggregated variable" % disaggregated_var.name) 1096 # There are two possibilities: It is declared on a Disjunct 1097 # transformation Block, or it is declared on the parent of a Disjunct 1098 # transformation block (if it is a single variable for multiple 1099 # Disjuncts the original doesn't appear in) 1100 transBlock = disaggregated_var.parent_block() 1101 if not hasattr(transBlock, '_disaggregatedVarMap'): 1102 try: 1103 transBlock = transBlock.parent_block().parent_block() 1104 except: 1105 logger.error(msg) 1106 raise 1107 try: 1108 return transBlock._disaggregatedVarMap['srcVar'][disaggregated_var] 1109 except: 1110 logger.error(msg) 1111 raise 1112 1113 # retrieves the disaggregation constraint for original_var resulting from 1114 # transforming disjunction 1115 def get_disaggregation_constraint(self, original_var, disjunction): 1116 """ 1117 Returns the disaggregation (re-aggregation?) constraint 1118 (which links the disaggregated variables to their original) 1119 corresponding to original_var and the transformation of disjunction. 1120 1121 Parameters 1122 ---------- 1123 original_var: a Var which was disaggregated in the transformation 1124 of Disjunction disjunction 1125 disjunction: a transformed Disjunction containing original_var 1126 """ 1127 for disjunct in disjunction.disjuncts: 1128 transBlock = disjunct._transformation_block 1129 if transBlock is not None: 1130 break 1131 if transBlock is None: 1132 raise GDP_Error("Disjunction '%s' has not been properly transformed:" 1133 " None of its disjuncts are transformed." 1134 % disjunction.name) 1135 1136 try: 1137 return transBlock().parent_block()._disaggregationConstraintMap[ 1138 original_var][disjunction] 1139 except: 1140 logger.error("It doesn't appear that '%s' is a variable that was " 1141 "disaggregated by Disjunction '%s'" % 1142 (original_var.name, disjunction.name)) 1143 raise 1144 1145 def get_var_bounds_constraint(self, v): 1146 """ 1147 Returns the IndexedConstraint which sets a disaggregated 1148 variable to be within its bounds when its Disjunct is active and to 1149 be 0 otherwise. (It is always an IndexedConstraint because each 1150 bound becomes a separate constraint.) 1151 1152 Parameters 1153 ---------- 1154 v: a Var which was created by the hull transformation as a 1155 disaggregated variable (and so appears on a transformation 1156 block of some Disjunct) 1157 """ 1158 msg = ("Either '%s' is not a disaggregated variable, or " 1159 "the disjunction that disaggregates it has not " 1160 "been properly transformed." % v.name) 1161 # This can only go well if v is a disaggregated var 1162 transBlock = v.parent_block() 1163 if not hasattr(transBlock, '_bigMConstraintMap'): 1164 try: 1165 transBlock = transBlock.parent_block().parent_block() 1166 except: 1167 logger.error(msg) 1168 raise 1169 try: 1170 return transBlock._bigMConstraintMap[v] 1171 except: 1172 logger.error(msg) 1173 raise 1174 1175 1176@TransformationFactory.register( 1177 'gdp.chull', 1178 doc="[DEPRECATED] please use 'gdp.hull' to get the Hull transformation.") 1179@deprecated("The 'gdp.chull' name is deprecated. " 1180 "Please use the more apt 'gdp.hull' instead.", 1181 logger='pyomo.gdp', 1182 version="5.7") 1183class _Deprecated_Name_Hull(Hull_Reformulation): 1184 def __init__(self): 1185 super(_Deprecated_Name_Hull, self).__init__() 1186