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 11from pyomo.core import (Var, Block, Constraint, Param, Set, SetOf, Suffix, 12 Expression, Objective, SortComponents, value, 13 ConstraintList) 14from pyomo.core.base import TransformationFactory, _VarData 15from pyomo.core.plugins.transform.hierarchy import Transformation 16from pyomo.common.config import ConfigBlock, ConfigValue, NonNegativeFloat 17from pyomo.common.modeling import unique_component_name 18from pyomo.repn.standard_repn import generate_standard_repn 19from pyomo.common.collections import ComponentMap, ComponentSet 20from pyomo.opt import TerminationCondition 21 22import logging 23 24logger = logging.getLogger('pyomo.contrib.fme') 25NAME_BUFFER = {} 26 27def _check_var_bounds_filter(constraint): 28 """Check if the constraint is already implied by the variable bounds""" 29 # this is one of our constraints, so we know that it is >=. 30 min_lhs = 0 31 for v, coef in constraint['map'].items(): 32 if coef > 0: 33 if v.lb is None: 34 return True # we don't have var bounds with which to imply the 35 # constraint... 36 min_lhs += coef*v.lb 37 elif coef < 0: 38 if v.ub is None: 39 return True # we don't have var bounds with which to imply the 40 # constraint... 41 min_lhs += coef*v.ub 42 # we do need value here since we didn't control v.lb and v.ub above. 43 if value(min_lhs) >= constraint['lower']: 44 return False # constraint implied by var bounds 45 return True 46 47def vars_to_eliminate_list(x): 48 if isinstance(x, (Var, _VarData)): 49 if not x.is_indexed(): 50 return ComponentSet([x]) 51 ans = ComponentSet() 52 for j in x.index_set(): 53 ans.add(x[j]) 54 return ans 55 elif hasattr(x, '__iter__'): 56 ans = ComponentSet() 57 for i in x: 58 ans.update(vars_to_eliminate_list(i)) 59 return ans 60 else: 61 raise ValueError( 62 "Expected Var or list of Vars." 63 "\n\tReceived %s" % type(x)) 64 65def gcd(a,b): 66 while b != 0: 67 a, b = b, a % b 68 return abs(a) 69 70def lcm(ints): 71 a = ints[0] 72 for b in ints[1:]: 73 a = abs(a*b) // gcd(a,b) 74 return a 75 76@TransformationFactory.register('contrib.fourier_motzkin_elimination', 77 doc="Project out specified (continuous) " 78 "variables from a linear model.") 79class Fourier_Motzkin_Elimination_Transformation(Transformation): 80 """Project out specified variables from a linear model. 81 82 This transformation requires the following keyword argument: 83 vars_to_eliminate: A user-specified list of continuous variables to 84 project out of the model 85 86 The transformation will deactivate the original constraints of the model 87 and create a new block named "_pyomo_contrib_fme_transformation" with the 88 projected constraints. Note that this transformation will flatten the 89 structure of the original model since there is no obvious mapping between 90 the original model and the transformed one. 91 92 """ 93 94 CONFIG = ConfigBlock("contrib.fourier_motzkin_elimination") 95 CONFIG.declare('vars_to_eliminate', ConfigValue( 96 default=None, 97 domain=vars_to_eliminate_list, 98 description="Continuous variable or list of continuous variables to " 99 "project out of the model", 100 doc=""" 101 This specifies the list of variables to project out of the model. 102 Note that these variables must all be continuous and the model must be 103 linear.""" 104 )) 105 CONFIG.declare('constraint_filtering_callback', ConfigValue( 106 default=_check_var_bounds_filter, 107 description="A callback that determines whether or not new " 108 "constraints generated by Fourier-Motzkin elimination are added " 109 "to the model", 110 doc=""" 111 Specify None in order for no constraint filtering to occur during the 112 transformation. 113 114 Specify a function that accepts a constraint (represented in the >= 115 dictionary form used in this transformation) and returns a Boolean 116 indicating whether or not to add it to the model. 117 """ 118 )) 119 CONFIG.declare('do_integer_arithmetic', ConfigValue( 120 default=False, 121 domain=bool, 122 description="A Boolean flag to decide whether Fourier-Motzkin " 123 "elimination will be performed with only integer arithmetic.", 124 doc=""" 125 If True, only integer arithmetic will be performed during Fourier- 126 Motzkin elimination. This should result in no numerical error. 127 If True and there is non-integer data in the constraints being 128 projected, an error will be raised. 129 130 If False, the algorithm will not check whether data is integer, and will 131 perform division operations. Use this setting when not all data is 132 integer, or when you are willing to sacrifice some numeric accuracy. 133 """ 134 )) 135 CONFIG.declare('verbose', ConfigValue( 136 default=False, 137 domain=bool, 138 description="A Boolean flag to enable verbose output.", 139 doc=""" 140 If True, logs the steps of the projection. 141 """ 142 )) 143 CONFIG.declare('zero_tolerance', ConfigValue( 144 default=0, 145 domain=NonNegativeFloat, 146 description="Absolute tolerance at which a float will be considered 0.", 147 doc=""" 148 Whenever fourier-motzkin elimination is used with non-integer data, 149 there is a chance of numeric trouble, the most obvious of which is 150 that 'eliminated' variables will remain in the constraints with very 151 small coefficients. Set this tolerance so that floating points smaller 152 than this will be treated as 0 (and reported that way in the final 153 constraints). 154 """ 155 )) 156 CONFIG.declare('integer_tolerance', ConfigValue( 157 default=0, 158 domain=NonNegativeFloat, 159 description="Absolute tolerance at which a float will be considered " 160 "(and cast to) an integer, when do_integer_arithmetic is True", 161 doc=""" 162 Tolerance at which a number x will be considered an integer, when we 163 are performing fourier-motzkin elimination with only integer_arithmetic. 164 That is, x will be cast to an integer if 165 abs(int(x) - x) <= integer_tolerance. 166 """ 167 )) 168 CONFIG.declare('projected_constraints_name', ConfigValue( 169 default=None, 170 domain=str, 171 description="Optional name for the ConstraintList containing the " 172 "projected constraints. Must be a unique name with respect to the " 173 "instance.", 174 doc=""" 175 Optional name for the ConstraintList containing the projected 176 constraints. If not specified, the constraints will be stored on a 177 private block created by the transformation, so if you want access 178 to them after the transformation, use this argument. 179 180 Must be a string which is a unique component name with respect to the 181 Block on which the transformation is called. 182 """ 183 )) 184 185 def __init__(self): 186 """Initialize transformation object""" 187 super(Fourier_Motzkin_Elimination_Transformation, self).__init__() 188 189 def _apply_to(self, instance, **kwds): 190 log_level = logger.level 191 try: 192 assert not NAME_BUFFER 193 config = self.CONFIG(kwds.pop('options', {})) 194 config.set_value(kwds) 195 # lower logging values emit more 196 if config.verbose and log_level > logging.INFO: 197 logger.setLevel(logging.INFO) 198 self.verbose = True 199 # if the user used the logger to ask for info level messages 200 elif log_level <= logging.INFO: 201 self.verbose = True 202 else: 203 self.verbose = False 204 self._apply_to_impl(instance, config) 205 finally: 206 # clear the global name buffer 207 NAME_BUFFER.clear() 208 # restore logging level 209 logger.setLevel(log_level) 210 211 def _apply_to_impl(self, instance, config): 212 vars_to_eliminate = config.vars_to_eliminate 213 self.constraint_filter = config.constraint_filtering_callback 214 self.do_integer_arithmetic = config.do_integer_arithmetic 215 self.integer_tolerance = config.integer_tolerance 216 self.zero_tolerance = config.zero_tolerance 217 if vars_to_eliminate is None: 218 raise RuntimeError("The Fourier-Motzkin Elimination transformation " 219 "requires the argument vars_to_eliminate, a " 220 "list of Vars to be projected out of the model.") 221 222 # make transformation block 223 transBlockName = unique_component_name( 224 instance, 225 '_pyomo_contrib_fme_transformation') 226 transBlock = Block() 227 instance.add_component(transBlockName, transBlock) 228 nm = config.projected_constraints_name 229 if nm is None: 230 projected_constraints = transBlock.projected_constraints = \ 231 ConstraintList() 232 else: 233 # check that this component doesn't already exist 234 if instance.component(nm) is not None: 235 raise RuntimeError("projected_constraints_name was specified " 236 "as '%s', but this is already a component " 237 "on the instance! Please specify a unique " 238 "name." % nm) 239 projected_constraints = ConstraintList() 240 instance.add_component(nm, projected_constraints) 241 242 # collect all of the constraints 243 # NOTE that we are ignoring deactivated constraints 244 constraints = [] 245 ctypes_not_to_transform = set((Block, Param, Objective, Set, SetOf, 246 Expression, Suffix, Var)) 247 for obj in instance.component_data_objects( 248 descend_into=Block, 249 sort=SortComponents.deterministic, 250 active=True): 251 if obj.ctype in ctypes_not_to_transform: 252 continue 253 elif obj.ctype is Constraint: 254 cons_list = self._process_constraint(obj) 255 constraints.extend(cons_list) 256 obj.deactivate() # the truth will be on our transformation block 257 else: 258 raise RuntimeError( 259 "Found active component %s of type %s. The " 260 "Fourier-Motzkin Elimination transformation can only " 261 "handle purely algebraic models. That is, only " 262 "Sets, Params, Vars, Constraints, Expressions, Blocks, " 263 "and Objectives may be active on the model." % (obj.name, 264 obj.ctype)) 265 266 for obj in vars_to_eliminate: 267 if obj.lb is not None: 268 constraints.append({'body': generate_standard_repn(obj), 269 'lower': value(obj.lb), 270 'map': ComponentMap([(obj, 1)])}) 271 if obj.ub is not None: 272 constraints.append({'body': generate_standard_repn(-obj), 273 'lower': -value(obj.ub), 274 'map': ComponentMap([(obj, -1)])}) 275 276 new_constraints = self._fourier_motzkin_elimination( constraints, 277 vars_to_eliminate) 278 279 # put the new constraints on the transformation block 280 for cons in new_constraints: 281 if self.constraint_filter is not None: 282 try: 283 keep = self.constraint_filter(cons) 284 except: 285 logger.error("Problem calling constraint filter callback " 286 "on constraint with right-hand side %s and " 287 "body:\n%s" % (cons['lower'], 288 cons['body'].to_expression())) 289 raise 290 if not keep: 291 continue 292 lhs = cons['body'].to_expression(sort=True) 293 lower = cons['lower'] 294 assert type(lower) is int or type(lower) is float 295 if type(lhs >= lower) is bool: 296 if lhs >= lower: 297 continue 298 else: 299 # This would actually make a lot of sense in this case... 300 #projected_constraints.add(Constraint.Infeasible) 301 raise RuntimeError("Fourier-Motzkin found the model is " 302 "infeasible!") 303 else: 304 projected_constraints.add(lhs >= lower) 305 306 def _process_constraint(self, constraint): 307 """Transforms a pyomo Constraint object into a list of dictionaries 308 representing only >= constraints. That is, if the constraint has both an 309 ub and a lb, it is transformed into two constraints. Otherwise it is 310 flipped if it is <=. Each dictionary contains the keys 'lower', 311 and 'body' where, after the process, 'lower' will be a constant, and 312 'body' will be the standard repn of the body. (The constant will be 313 moved to the RHS and we know that the upper bound is None after this). 314 """ 315 body = constraint.body 316 std_repn = generate_standard_repn(body) 317 # make sure that we store the lower bound's value so that we need not 318 # worry again during the transformation 319 cons_dict = {'lower': value(constraint.lower), 320 'body': std_repn 321 } 322 upper = value(constraint.upper) 323 constraints_to_add = [cons_dict] 324 if upper is not None: 325 # if it has both bounds 326 if cons_dict['lower'] is not None: 327 # copy the constraint and flip 328 leq_side = {'lower': -upper, 329 'body': generate_standard_repn(-1.0*body)} 330 self._move_constant_and_add_map(leq_side) 331 constraints_to_add.append(leq_side) 332 333 # If it has only an upper bound, we just need to flip it 334 else: 335 # just flip the constraint 336 cons_dict['lower'] = -upper 337 cons_dict['body'] = generate_standard_repn(-1.0*body) 338 self._move_constant_and_add_map(cons_dict) 339 340 return constraints_to_add 341 342 def _move_constant_and_add_map(self, cons_dict): 343 """Takes constraint in dicionary form already in >= form, 344 and moves the constant to the RHS 345 """ 346 body = cons_dict['body'] 347 constant = value(body.constant) 348 cons_dict['lower'] -= constant 349 body.constant = 0 350 351 # store a map of vars to coefficients. We can't use this in place of 352 # standard repn because determinism, but this will save a lot of linear 353 # time searches later. Note also that we will take the value of the 354 # coeficient here so that we never have to worry about it again during 355 # the transformation. 356 cons_dict['map'] = ComponentMap(zip(body.linear_vars, 357 [value(coef) for coef in 358 body.linear_coefs])) 359 360 def _fourier_motzkin_elimination(self, constraints, vars_to_eliminate): 361 """Performs FME on the constraint list in the argument 362 (which is assumed to be all >= constraints and stored in the 363 dictionary representation), projecting out each of the variables in 364 vars_to_eliminate""" 365 366 # We only need to eliminate variables that actually appear in 367 # this set of constraints... Revise our list. 368 vars_that_appear = [] 369 vars_that_appear_set = ComponentSet() 370 for cons in constraints: 371 std_repn = cons['body'] 372 if not std_repn.is_linear(): 373 # as long as none of vars_that_appear are in the nonlinear part, 374 # we are actually okay. 375 nonlinear_vars = ComponentSet(v for two_tuple in 376 std_repn.quadratic_vars for 377 v in two_tuple) 378 nonlinear_vars.update(v for v in std_repn.nonlinear_vars) 379 for var in nonlinear_vars: 380 if var in vars_to_eliminate: 381 raise RuntimeError("Variable %s appears in a nonlinear " 382 "constraint. The Fourier-Motzkin " 383 "Elimination transformation can only " 384 "be used to eliminate variables " 385 "which only appear linearly." % 386 var.name) 387 for var in std_repn.linear_vars: 388 if var in vars_to_eliminate: 389 if not var in vars_that_appear_set: 390 vars_that_appear.append(var) 391 vars_that_appear_set.add(var) 392 393 # we actually begin the recursion here 394 total = len(vars_that_appear) 395 iteration = 1 396 while vars_that_appear: 397 # first var we will project out 398 the_var = vars_that_appear.pop() 399 logger.warning("Projecting out var %s of %s" % (iteration, total)) 400 if self.verbose: 401 logger.info("Projecting out %s" % 402 the_var.getname(fully_qualified=True, 403 name_buffer=NAME_BUFFER)) 404 logger.info("New constraints are:") 405 406 # we are 'reorganizing' the constraints, we sort based on the sign 407 # of the coefficient of the_var: This tells us whether we have 408 # the_var <= other stuff or vice versa. 409 leq_list = [] 410 geq_list = [] 411 waiting_list = [] 412 413 coefs = [] 414 for cons in constraints: 415 leaving_var_coef = cons['map'].get(the_var) 416 if leaving_var_coef is None or leaving_var_coef == 0: 417 waiting_list.append(cons) 418 if self.verbose: 419 logger.info("\t%s <= %s" 420 % (cons['lower'], 421 cons['body'].to_expression())) 422 continue 423 424 # we know the constraint is a >= constraint, using that 425 # assumption below. 426 # NOTE: neither of the scalar multiplications below flip the 427 # constraint. So we are sure to have only geq constraints 428 # forever, which is exactly what we want. 429 if not self.do_integer_arithmetic: 430 if leaving_var_coef < 0: 431 leq_list.append( 432 self._nonneg_scalar_multiply_linear_constraint( 433 cons, -1.0/leaving_var_coef)) 434 else: 435 geq_list.append( 436 self._nonneg_scalar_multiply_linear_constraint( 437 cons, 1.0/leaving_var_coef)) 438 else: 439 coefs.append(self._as_integer( 440 leaving_var_coef, 441 self._get_noninteger_coef_error_message, 442 (the_var.name, leaving_var_coef) 443 )) 444 if self.do_integer_arithmetic and len(coefs) > 0: 445 least_common_mult = lcm(coefs) 446 for cons in constraints: 447 leaving_var_coef = cons['map'].get(the_var) 448 if leaving_var_coef is None or leaving_var_coef == 0: 449 continue 450 to_lcm = least_common_mult // abs(int(leaving_var_coef)) 451 if leaving_var_coef < 0: 452 leq_list.append( 453 self._nonneg_scalar_multiply_linear_constraint( 454 cons, to_lcm)) 455 else: 456 geq_list.append( 457 self._nonneg_scalar_multiply_linear_constraint( 458 cons, to_lcm)) 459 460 constraints = waiting_list 461 for leq in leq_list: 462 for geq in geq_list: 463 constraints.append( self._add_linear_constraints( leq, geq)) 464 if self.verbose: 465 cons = constraints[len(constraints)-1] 466 logger.info("\t%s <= %s" % 467 (cons['lower'], 468 cons['body'].to_expression())) 469 470 iteration += 1 471 472 return constraints 473 474 def _get_noninteger_coef_error_message(self, varname, coef): 475 return ("The do_integer_arithmetic flag was " 476 "set to True, but the coefficient of " 477 "%s is non-integer within the specified " 478 "tolerance, with value %s. \n" 479 "Please set do_integer_arithmetic=" 480 "False, increase integer_tolerance, " 481 "or make your data integer." % (varname, coef)) 482 483 def _as_integer(self, x, error_message, error_args): 484 if abs(int(x) - x) <= self.integer_tolerance: 485 return int(round(x)) 486 raise ValueError(error_message if error_args is None 487 else error_message(*error_args)) 488 489 def _multiply(self, scalar, coef, error_message, error_args): 490 if self.do_integer_arithmetic: 491 assert type(scalar) is int 492 return scalar * self._as_integer(coef, error_message, error_args) 493 elif abs(scalar*coef) > self.zero_tolerance: 494 return scalar*coef 495 else: 496 return 0 497 498 def _add(self, a, b, error_message, error_args): 499 if self.do_integer_arithmetic: 500 return self._as_integer(a, error_message, error_args) \ 501 + self._as_integer(b, error_message, error_args) 502 elif abs(a + b) > self.zero_tolerance: 503 return a + b 504 else: 505 return 0 506 507 def _nonneg_scalar_multiply_linear_constraint_error_msg(self, cons, coef): 508 return ( 509 "The do_integer_arithmetic flag was set to True, but the " 510 "lower bound of %s is non-integer within the specified " 511 "tolerance, with value %s. \n" 512 "Please set do_integer_arithmetic=False, increase " 513 "integer_tolerance, or make your data integer." % 514 (cons['body'].to_expression() >= cons['lower'], coef) 515 ) 516 517 def _nonneg_scalar_multiply_linear_constraint(self, cons, scalar): 518 """Multiplies all coefficients and the RHS of a >= constraint by scalar. 519 There is no logic for flipping the equality, so this is just the 520 special case with a nonnegative scalar, which is all we need. 521 522 If self.do_integer_arithmetic is True, this assumes that scalar is an 523 int. It also will throw an error if any data is non-integer (within 524 tolerance) 525 """ 526 body = cons['body'] 527 new_coefs = [] 528 for i, coef in enumerate(body.linear_coefs): 529 v = body.linear_vars[i] 530 new_coefs.append(self._multiply( 531 scalar, coef, self._get_noninteger_coef_error_message, 532 (v.name, coef) 533 )) 534 # update the map 535 cons['map'][v] = new_coefs[i] 536 body.linear_coefs = new_coefs 537 538 body.quadratic_coefs = [scalar*coef for coef in body.quadratic_coefs] 539 body.nonlinear_expr = scalar*body.nonlinear_expr if \ 540 body.nonlinear_expr is not None else None 541 542 # assume scalar >= 0 and constraint only has lower bound 543 lb = cons['lower'] 544 if lb is not None: 545 cons['lower'] = self._multiply( 546 scalar, lb, 547 self._nonneg_scalar_multiply_linear_constraint_error_msg, 548 (cons, coef) 549 ) 550 return cons 551 552 def _add_linear_constraints_error_msg(self, cons1, cons2): 553 return ( 554 "The do_integer_arithmetic flag was set to True, but while " 555 "adding %s and %s, encountered a coefficient that is " 556 "non-integer within the specified tolerance\n" 557 "Please set do_integer_arithmetic=False, increase " 558 "integer_tolerance, or make your data integer." % 559 (cons1['body'].to_expression() >= cons1['lower'], 560 cons2['body'].to_expression() >= cons2['lower']) 561 ) 562 563 def _add_linear_constraints(self, cons1, cons2): 564 """Adds two >= constraints 565 566 Because this is always called after 567 _nonneg_scalar_multiply_linear_constraint, though it is implemented 568 more generally. 569 """ 570 ans = {'lower': None, 'body': None, 'map': ComponentMap()} 571 cons1_body = cons1['body'] 572 cons2_body = cons2['body'] 573 574 # Need this to be both deterministic and to account for the fact that 575 # Vars aren't hashable. 576 all_vars = list(cons1_body.linear_vars) 577 seen = ComponentSet(all_vars) 578 for v in cons2_body.linear_vars: 579 if v not in seen: 580 all_vars.append(v) 581 582 expr = 0 583 for var in all_vars: 584 coef = self._add( 585 cons1['map'].get(var, 0), cons2['map'].get(var, 0), 586 self._add_linear_constraints_error_msg, (cons1, cons2)) 587 ans['map'][var] = coef 588 expr += coef*var 589 590 # deal with nonlinear stuff if there is any 591 for cons in [cons1_body, cons2_body]: 592 if cons.nonlinear_expr is not None: 593 expr += cons.nonlinear_expr 594 expr += sum(coef*v1*v2 for (coef, (v1, v2)) in 595 zip(cons.quadratic_coefs, cons.quadratic_vars)) 596 597 ans['body'] = generate_standard_repn(expr) 598 599 # upper is None and lower exists, so this gets the constant 600 ans['lower'] = self._add( 601 cons1['lower'], cons2['lower'], 602 self._add_linear_constraints_error_msg, (cons1, cons2)) 603 604 return ans 605 606 def post_process_fme_constraints(self, m, solver_factory, 607 projected_constraints=None, tolerance=0): 608 """Function that solves a sequence of LPs problems to check if 609 constraints are implied by each other. Deletes any that are. 610 611 Parameters 612 ---------------- 613 m: A model, already transformed with FME. Note that if constraints 614 have been added, activated, or deactivated, we will check for 615 redundancy against the whole active part of the model. If you call 616 this straight after FME, you are only checking within the projected 617 constraints, but otherwise it is up to the user. 618 solver_factory: A SolverFactory object (constructed with a solver 619 which can solve the continuous relaxation of the 620 active constraints on the model. That is, if you 621 had nonlinear constraints unrelated to the variables 622 being projected, you need to either deactivate them or 623 provide a solver which will do the right thing.) 624 projected_constraints: The ConstraintList of projected constraints. 625 Default is None, in which case we assume that 626 the FME transformation was called without 627 specifying their name, so will look for them on 628 the private transformation block. 629 tolerance: Tolerance at which we decide a constraint is implied by the 630 others. Default is 0, meaning we remove the constraint if 631 the LP solve finds the constraint can be tight but not 632 violated. Setting this to a small positive value would 633 remove constraints more conservatively. Setting it to a 634 negative value would result in a relaxed problem. 635 """ 636 if projected_constraints is None: 637 # make sure m looks like what we expect 638 if not hasattr(m, "_pyomo_contrib_fme_transformation"): 639 raise RuntimeError("It looks like model %s has not been " 640 "transformed with the " 641 "fourier_motzkin_elimination transformation!" 642 % m.name) 643 transBlock = m._pyomo_contrib_fme_transformation 644 if not hasattr(transBlock, 'projected_constraints'): 645 raise RuntimeError("It looks the projected constraints " 646 "were manually named when the FME " 647 "transformation was called on %s. " 648 "If this is so, specify the ConstraintList " 649 "of projected constraints with the " 650 "'projected_constraints' argument." % m.name) 651 projected_constraints = transBlock.projected_constraints 652 653 # relax integrality so that we can do this with LP solves. 654 TransformationFactory('core.relax_integer_vars').apply_to( 655 m, transform_deactivated_blocks=True) 656 # deactivate any active objectives on the model, and save what we did so 657 # we can undo it after. 658 active_objs = [] 659 for obj in m.component_data_objects(Objective, descend_into=True): 660 if obj.active: 661 active_objs.append(obj) 662 obj.deactivate() 663 # add placeholder for our own objective 664 obj_name = unique_component_name(m, '_fme_post_process_obj') 665 obj = Objective(expr=0) 666 m.add_component(obj_name, obj) 667 for i in projected_constraints: 668 # If someone wants us to ignore it and leave it in the model, we 669 # can. 670 if not projected_constraints[i].active: 671 continue 672 # deactivate the constraint 673 projected_constraints[i].deactivate() 674 m.del_component(obj) 675 # make objective to maximize its infeasibility 676 obj = Objective(expr=projected_constraints[i].body - \ 677 projected_constraints[i].lower) 678 m.add_component(obj_name, obj) 679 results = solver_factory.solve(m) 680 if results.solver.termination_condition == \ 681 TerminationCondition.unbounded: 682 obj_val = -float('inf') 683 elif results.solver.termination_condition != \ 684 TerminationCondition.optimal: 685 raise RuntimeError("Unsuccessful subproblem solve when checking" 686 "constraint %s.\n\t" 687 "Termination Condition: %s" % 688 (projected_constraints[i].name, 689 results.solver.termination_condition)) 690 else: 691 obj_val = value(obj) 692 # if we couldn't make it infeasible, it's useless 693 if obj_val >= tolerance: 694 m.del_component(projected_constraints[i]) 695 del projected_constraints[i] 696 else: 697 projected_constraints[i].activate() 698 699 # clean up 700 m.del_component(obj) 701 for obj in active_objs: 702 obj.activate() 703 # undo relax integrality 704 TransformationFactory('core.relax_integer_vars').apply_to(m, undo=True) 705