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""" 12Cutting plane-based GDP reformulation. 13 14Implements a general cutting plane-based reformulation for linear and 15convex GDPs. 16""" 17from __future__ import division 18 19from pyomo.common.config import (ConfigBlock, ConfigValue, 20 NonNegativeFloat, PositiveInt, In) 21from pyomo.common.modeling import unique_component_name 22from pyomo.core import ( Any, Block, Constraint, Objective, Param, Var, 23 SortComponents, Transformation, TransformationFactory, 24 value, NonNegativeIntegers, Reals, NonNegativeReals, 25 Suffix, ComponentMap ) 26from pyomo.core.expr import differentiate 27from pyomo.common.collections import ComponentSet 28from pyomo.opt import SolverFactory 29from pyomo.repn import generate_standard_repn 30 31from pyomo.gdp import Disjunct, Disjunction, GDP_Error 32from pyomo.gdp.util import ( verify_successful_solve, NORMAL, 33 clone_without_expression_components ) 34 35from pyomo.contrib.fme.fourier_motzkin_elimination import \ 36 Fourier_Motzkin_Elimination_Transformation 37 38import logging 39 40logger = logging.getLogger('pyomo.gdp.cuttingplane') 41 42NAME_BUFFER = {} 43 44def do_not_tighten(m): 45 return m 46 47def _get_constraint_exprs(constraints, hull_to_bigm_map): 48 """Returns a list of expressions which are constrain.expr translated 49 into the bigm space, for each constraint in constraints. 50 """ 51 cuts = [] 52 for cons in constraints: 53 cuts.append(clone_without_expression_components( 54 cons.expr, substitute=hull_to_bigm_map)) 55 return cuts 56 57def _constraint_tight(model, constraint, TOL): 58 """ 59 Returns a list [a,b] where a is -1 if the lower bound is tight or 60 slightly violated, b is 1 if the upper bound is tight of slightly 61 violated, and [a,b]=[-1,1] if we have an exactly satisfied (or 62 slightly violated) equality. 63 """ 64 val = value(constraint.body) 65 ans = [0, 0] 66 if constraint.lower is not None: 67 if val - value(constraint.lower) <= TOL: 68 # tight or in violation of LB 69 ans[0] -= 1 70 71 if constraint.upper is not None: 72 if value(constraint.upper) - val <= TOL: 73 # tight or in violation of UB 74 ans[1] += 1 75 76 return ans 77 78def _get_linear_approximation_expr(normal_vec, point): 79 """Returns constraint linearly approximating constraint normal to normal_vec 80 at point""" 81 body = 0 82 for coef, v in zip(point, normal_vec): 83 body -= coef*v 84 return body >= -sum(normal_vec[idx]*v.value for (idx, v) in 85 enumerate(point)) 86 87def _precompute_potentially_useful_constraints(transBlock_rHull, 88 disaggregated_vars): 89 instance_rHull = transBlock_rHull.model() 90 constraints = transBlock_rHull.constraints_for_FME = [] 91 for constraint in instance_rHull.component_data_objects( 92 Constraint, 93 active=True, 94 descend_into=Block, 95 sort=SortComponents.deterministic): 96 # we don't care about anything that does not involve at least one 97 # disaggregated variable. 98 repn = generate_standard_repn(constraint.body) 99 for v in repn.linear_vars + repn.quadratic_vars + repn.nonlinear_vars: 100 # ESJ: This is why disaggregated_vars is a ComponentSet 101 if v in disaggregated_vars: 102 constraints.append(constraint) 103 break 104 105def create_cuts_fme(transBlock_rHull, var_info, hull_to_bigm_map, 106 rBigM_linear_constraints, rHull_vars, disaggregated_vars, 107 norm, cut_threshold, zero_tolerance, integer_arithmetic, 108 constraint_tolerance): 109 """Returns a cut which removes x* from the relaxed bigm feasible region. 110 111 Finds all the constraints which are tight at xhat (assumed to be the 112 solution currently in instance_rHull), and calculates a composite normal 113 vector by summing the vectors normal to each of these constraints. Then 114 Fourier-Motzkin elimination is used to project the disaggregated variables 115 out of the polyhedron formed by the composite normal and the collection 116 of tight constraints. This results in multiple cuts, of which we select 117 one that cuts of x* by the greatest margin, as long as that margin is 118 more than cut_threshold. If no cut satisfies the margin specified by 119 cut_threshold, we return None. 120 121 Parameters 122 ----------- 123 transBlock_rHull: transformation blcok on relaxed hull instance 124 var_info: List of tuples (rBigM_var, rHull_var, xstar_param) 125 hull_to_bigm_map: For expression substition, maps id(hull_var) to 126 coresponding bigm var 127 rBigM_linear_constraints: list of linear constraints in relaxed bigM 128 rHull_vars: list of all variables in relaxed hull 129 disaggregated_vars: ComponentSet of disaggregated variables in hull 130 reformulation 131 norm: norm used in the separation problem 132 cut_threshold: Amount x* needs to be infeasible in generated cut in order 133 to consider the cut for addition to the bigM model. 134 zero_tolerance: Tolerance at which a float will be treated as 0 during 135 Fourier-Motzkin elimination 136 integer_arithmetic: boolean, whether or not to require Fourier-Motzkin 137 Elimination does integer arithmetic. Only possible 138 when all data is integer. 139 constraint_tolerance: Tolerance at which we will consider a constraint 140 tight. 141 """ 142 instance_rHull = transBlock_rHull.model() 143 # In the first iteration, we will compute a list of constraints that could 144 # ever be interesting: Everything that involves at least one disaggregated 145 # variable. 146 if transBlock_rHull.component("constraints_for_FME") is None: 147 _precompute_potentially_useful_constraints( transBlock_rHull, 148 disaggregated_vars) 149 150 tight_constraints = Block() 151 conslist = tight_constraints.constraints = Constraint( 152 NonNegativeIntegers) 153 conslist.construct() 154 something_interesting = False 155 for constraint in transBlock_rHull.constraints_for_FME: 156 multipliers = _constraint_tight(instance_rHull, constraint, 157 constraint_tolerance) 158 for multiplier in multipliers: 159 if multiplier: 160 something_interesting = True 161 f = constraint.body 162 firstDerivs = differentiate(f, wrt_list=rHull_vars) 163 normal_vec = [multiplier*value(_) for _ in firstDerivs] 164 # check if constraint is linear 165 if f.polynomial_degree() == 1: 166 conslist[len(conslist)] = constraint.expr 167 else: 168 # we will use the linear approximation of this constraint at 169 # x_hat 170 conslist[len(conslist)] = _get_linear_approximation_expr( 171 normal_vec, rHull_vars) 172 173 # NOTE: we now have all the tight Constraints (in the pyomo sense of the 174 # word "Constraint"), but we are missing some variable bounds. The ones for 175 # the disaggregated variables will be added by FME 176 177 # It is possible that the separation problem returned a point in the 178 # interior of the convex hull. It is also possible that the only active 179 # constraints do not involve the disaggregated variables. In these 180 # situations, there are not constraints from which to create a valid cut. 181 if not something_interesting: 182 return None 183 184 tight_constraints.construct() 185 logger.info("Calling FME transformation on %s constraints to eliminate" 186 " %s variables" % (len(tight_constraints.constraints), 187 len(disaggregated_vars))) 188 TransformationFactory('contrib.fourier_motzkin_elimination').\ 189 apply_to(tight_constraints, vars_to_eliminate=disaggregated_vars, 190 zero_tolerance=zero_tolerance, 191 do_integer_arithmetic=integer_arithmetic, 192 projected_constraints_name="fme_constraints") 193 fme_results = tight_constraints.fme_constraints 194 projected_constraints = [cons for i, cons in fme_results.items()] 195 196 # we created these constraints with the variables from rHull. We 197 # actually need constraints for BigM and rBigM now! 198 cuts = _get_constraint_exprs(projected_constraints, hull_to_bigm_map) 199 200 # We likely have some cuts that duplicate other constraints now. We will 201 # filter them to make sure that they do in fact cut off x*. If that's the 202 # case, we know they are not already in the BigM relaxation. Because they 203 # came from FME, they are very likely redundant, so we'll keep the best one 204 # we find 205 best = 0 206 best_cut = None 207 cuts_to_keep = [] 208 for i, cut in enumerate(cuts): 209 # x* is still in rBigM, so we can just remove this constraint if it 210 # is satisfied at x* 211 logger.info("FME: Post-processing cut %s" % cut) 212 if value(cut): 213 logger.info("FME:\t Doesn't cut off x*") 214 continue 215 # we have found a constraint which cuts of x* by some convincing amount 216 # and is not already in rBigM. 217 cuts_to_keep.append(i) 218 # We know cut is lb <= expr and that it's violated 219 assert len(cut.args) == 2 220 cut_off = value(cut.args[0]) - value(cut.args[1]) 221 if cut_off > cut_threshold and cut_off > best: 222 best = cut_off 223 best_cut = cut 224 logger.info("FME:\t New best cut: Cuts off x* by %s." % best) 225 226 # NOTE: this is not used right now, but it's not hard to imagine a world in 227 # which we would want to keep multiple cuts from FME, so leaving it in for 228 # now. 229 cuts = [cuts[i] for i in cuts_to_keep] 230 231 if best_cut is not None: 232 return [best_cut] 233 234 return None 235 236def create_cuts_normal_vector(transBlock_rHull, var_info, hull_to_bigm_map, 237 rBigM_linear_constraints, rHull_vars, 238 disaggregated_vars, norm, cut_threshold, 239 zero_tolerance, integer_arithmetic, 240 constraint_tolerance): 241 """Returns a cut which removes x* from the relaxed bigm feasible region. 242 243 Ignores all parameters except var_info and cut_threshold, and constructs 244 a cut at x_hat, the projection of the relaxed bigM solution x* onto the hull, 245 which is perpendicular to the vector from x* to x_hat. 246 247 Note that this method will often lead to numerical difficulties since both 248 x* and x_hat are solutions to optimization problems. To mitigate this, 249 use some method of backing off the cut to make it a bit more conservative. 250 251 Parameters 252 ----------- 253 transBlock_rHull: transformation blcok on relaxed hull instance. Ignored by 254 this callback. 255 var_info: List of tuples (rBigM_var, rHull_var, xstar_param) 256 hull_to_bigm_map: For expression substition, maps id(hull_var) to 257 coresponding bigm var. Ignored by this callback 258 rBigM_linear_constraints: list of linear constraints in relaxed bigM. 259 Ignored by this callback. 260 rHull_vars: list of all variables in relaxed hull. Ignored by this callback. 261 disaggregated_vars: ComponentSet of disaggregated variables in hull 262 reformulation. Ignored by this callback 263 norm: The norm used in the separation problem, will be used to calculate 264 the subgradient used to generate the cut 265 cut_threshold: Amount x* needs to be infeasible in generated cut in order 266 to consider the cut for addition to the bigM model. 267 zero_tolerance: Tolerance at which a float will be treated as 0 during 268 Fourier-Motzkin elimination. Ignored by this callback 269 integer_arithmetic: Ignored by this callback (specifies FME use integer 270 arithmetic) 271 constraint_tolerance: Ignored by this callback (specifies when constraints 272 are considered tight in FME) 273 """ 274 cutexpr = 0 275 if norm == 2: 276 for x_rbigm, x_hull, x_star in var_info: 277 cutexpr += (x_hull.value - x_star.value)*(x_rbigm - x_hull.value) 278 elif norm == float('inf'): 279 duals = transBlock_rHull.model().dual 280 if len(duals) == 0: 281 raise GDP_Error("No dual information in the separation problem! " 282 "To use the infinity norm and the " 283 "create_cuts_normal_vector method, you must use " 284 "a solver which provides dual information.") 285 i = 0 286 for x_rbigm, x_hull, x_star in var_info: 287 # ESJ: We wrote this so duals will be nonnegative 288 mu_plus = value(duals[transBlock_rHull.inf_norm_linearization[i]]) 289 mu_minus = value(duals[transBlock_rHull.inf_norm_linearization[i+1]]) 290 assert mu_plus >= 0 291 assert mu_minus >= 0 292 cutexpr += (mu_plus - mu_minus)*(x_rbigm - x_hull.value) 293 i += 2 294 295 # make sure we're cutting off x* by enough. 296 if value(cutexpr) < -cut_threshold: 297 return [cutexpr >= 0] 298 logger.warning("Generated cut did not remove relaxed BigM solution by more " 299 "than the specified threshold of %s. Stopping cut " 300 "generation." % cut_threshold) 301 return None 302 303def back_off_constraint_with_calculated_cut_violation(cut, transBlock_rHull, 304 bigm_to_hull_map, opt, 305 stream_solver, TOL): 306 """Calculates the maximum violation of cut subject to the relaxed hull 307 constraints. Increases this violation by TOL (to account for optimality 308 tolerance in solving the problem), and, if it finds that cut can be violated 309 up to this tolerance, makes it more conservative such that it no longer can. 310 311 Parameters 312 ---------- 313 cut: The cut to be made more conservative, a Constraint 314 transBlock_rHull: the relaxed hull model's transformation Block 315 bigm_to_hull_map: Dictionary mapping ids of bigM variables to the 316 corresponding variables on the relaxed hull instance 317 opt: SolverFactory object for solving the maximum violation problem 318 stream_solver: Whether or not to set tee=True while solving the maximum 319 violation problem. 320 TOL: An absolute tolerance to be added to the calculated cut violation, 321 to account for optimality tolerance in the maximum violation problem 322 solve. 323 """ 324 instance_rHull = transBlock_rHull.model() 325 logger.info("Post-processing cut: %s" % cut.expr) 326 # Take a constraint. We will solve a problem maximizing its violation 327 # subject to rHull. We will add some user-specified tolerance to that 328 # violation, and then add that much padding to it if it can be violated. 329 transBlock_rHull.separation_objective.deactivate() 330 331 transBlock_rHull.infeasibility_objective = Objective( 332 expr=clone_without_expression_components(cut.body, 333 substitute=bigm_to_hull_map)) 334 335 results = opt.solve(instance_rHull, tee=stream_solver, load_solutions=False) 336 if verify_successful_solve(results) is not NORMAL: 337 logger.warning("Problem to determine how much to " 338 "back off the new cut " 339 "did not solve normally. Leaving the constraint as is, " 340 "which could lead to numerical trouble%s" % (results,)) 341 # restore the objective 342 transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective) 343 transBlock_rHull.separation_objective.activate() 344 return 345 instance_rHull.solutions.load_from(results) 346 347 # we're minimizing, val is <= 0 348 val = value(transBlock_rHull.infeasibility_objective) - TOL 349 if val <= 0: 350 logger.info("\tBacking off cut by %s" % val) 351 cut._body += abs(val) 352 # else there is nothing to do: restore the objective 353 transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective) 354 transBlock_rHull.separation_objective.activate() 355 356def back_off_constraint_by_fixed_tolerance(cut, transBlock_rHull, 357 bigm_to_hull_map, opt, stream_solver, 358 TOL): 359 """Makes cut more conservative by absolute tolerance TOL 360 361 Parameters 362 ---------- 363 cut: the cut to be made more conservative, a Constraint 364 transBlock_rHull: the relaxed hull model's transformation Block. Ignored by 365 this callback 366 bigm_to_hull_map: Dictionary mapping ids of bigM variables to the 367 corresponding variables on the relaxed hull instance. 368 Ignored by this callback. 369 opt: SolverFactory object. Ignored by this callback 370 stream_solver: Whether or not to set tee=True while solving. Ignored by 371 this callback 372 TOL: An absolute tolerance to be added to make cut more conservative. 373 """ 374 cut._body += TOL 375 376@TransformationFactory.register('gdp.cuttingplane', 377 doc="Relaxes a linear disjunctive model by " 378 "adding cuts from convex hull to Big-M " 379 "reformulation.") 380class CuttingPlane_Transformation(Transformation): 381 """Relax convex disjunctive model by forming the bigm relaxation and then 382 iteratively adding cuts from the hull relaxation (or the hull relaxation 383 after some basic steps) in order to strengthen the formulation. 384 385 Note that gdp.cuttingplane is not a structural transformation: If variables 386 on the model are fixed, they will be treated as data, and unfixing them 387 after transformation will very likely result in an invalid model. 388 389 This transformation accepts the following keyword arguments: 390 391 Parameters 392 ---------- 393 solver : Solver name (as string) to use to solve relaxed BigM and separation 394 problems 395 solver_options : dictionary of options to pass to the solver 396 stream_solver : Whether or not to display solver output 397 verbose : Enable verbose output from cuttingplanes algorithm 398 cuts_name : Optional name for the IndexedConstraint containing the projected 399 cuts (must be a unique name with respect to the instance) 400 minimum_improvement_threshold : Stopping criterion based on improvement in 401 Big-M relaxation. This is the minimum 402 difference in relaxed BigM objective 403 values between consecutive iterations 404 separation_objective_threshold : Stopping criterion based on separation 405 objective. If separation objective is not 406 at least this large, cut generation will 407 terminate. 408 cut_filtering_threshold : Stopping criterion based on effectiveness of the 409 generated cut: This is the amount by which 410 a cut must be violated at the relaxed bigM 411 solution in order to be added to the bigM model 412 max_number_of_cuts : The maximum number of cuts to add to the big-M model 413 norm : norm to use in the objective of the separation problem 414 tighten_relaxation : callback to modify the GDP model before the hull 415 relaxation is taken (e.g. could be used to perform 416 basic steps) 417 create_cuts : callback to create cuts using the solved relaxed bigM and hull 418 problems 419 post_process_cut : callback to perform post-processing on created cuts 420 back_off_problem_tolerance : tolerance to use while post-processing 421 zero_tolerance : Tolerance at which a float will be considered 0 when 422 using Fourier-Motzkin elimination to create cuts. 423 do_integer_arithmetic : Whether or not to require Fourier-Motzkin elimination 424 to do integer arithmetic. Only possible when all 425 data is integer. 426 tight_constraint_tolerance : Tolerance at which a constraint is considered 427 tight for the Fourier-Motzkin cut generation 428 procedure 429 430 By default, the callbacks will be set such that the algorithm performed is 431 as presented in [1], but with an additional post-processing procedure to 432 reduce numerical error, which calculates the maximum violation of the cut 433 subject to the relaxed hull constraints, and then pads the constraint by 434 this violation plus an additional user-specified tolerance. 435 436 In addition, the create_cuts_fme function provides an (exponential time) 437 method of generating cuts which reduces numerical error (and can eliminate 438 it if all data is integer). It collects the hull constraints which are 439 tight at the solution of the separation problem. It creates a cut in the 440 extended space perpendicular to a composite normal vector created by 441 summing the directions normal to these constraints. It then performs 442 fourier-motzkin elimination on the collection of constraints and the cut 443 to project out the disaggregated variables. The resulting constraint which 444 is most violated by the relaxed bigM solution is then returned. 445 446 References 447 ---------- 448 [1] Sawaya, N. W., Grossmann, I. E. (2005). A cutting plane method for 449 solving linear generalized disjunctive programming problems. Computers 450 and Chemical Engineering, 29, 1891-1913 451 """ 452 453 CONFIG = ConfigBlock("gdp.cuttingplane") 454 CONFIG.declare('solver', ConfigValue( 455 default='ipopt', 456 domain=str, 457 description="""Solver to use for relaxed BigM problem and the separation 458 problem""", 459 doc=""" 460 This specifies the solver which will be used to solve LP relaxation 461 of the BigM problem and the separation problem. Note that this solver 462 must be able to handle a quadratic objective because of the separation 463 problem. 464 """ 465 )) 466 CONFIG.declare('minimum_improvement_threshold', ConfigValue( 467 default=0.01, 468 domain=NonNegativeFloat, 469 description="Threshold value for difference in relaxed bigM problem " 470 "objectives used to decide when to stop adding cuts", 471 doc=""" 472 If the difference between the objectives in two consecutive iterations is 473 less than this value, the algorithm terminates without adding the cut 474 generated in the last iteration. 475 """ 476 )) 477 CONFIG.declare('separation_objective_threshold', ConfigValue( 478 default=0.01, 479 domain=NonNegativeFloat, 480 description="Threshold value used to decide when to stop adding cuts: " 481 "If separation problem objective is not at least this quantity, cut " 482 "generation will terminate.", 483 doc=""" 484 If the separation problem objective (distance between relaxed bigM 485 solution and its projection onto the relaxed hull feasible region) 486 does not exceed this threshold, the algorithm will terminate. 487 """ 488 )) 489 CONFIG.declare('max_number_of_cuts', ConfigValue( 490 default=100, 491 domain=PositiveInt, 492 description="The maximum number of cuts to add before the algorithm " 493 "terminates.", 494 doc=""" 495 If the algorithm does not terminate due to another criterion first, 496 cut generation will stop after adding this many cuts. 497 """ 498 )) 499 CONFIG.declare('norm', ConfigValue( 500 default=2, 501 domain=In([2, float('inf')]), 502 description="Norm to use in the separation problem: 2, or " 503 "float('inf')", 504 doc=""" 505 Norm used to calculate distance in the objective of the separation 506 problem which finds the nearest point on the hull relaxation region 507 to the current solution of the relaxed bigm problem. 508 509 Supported norms are the Euclidean norm (specify 2) and the infinity 510 norm (specify float('inf')). Note that the first makes the separation 511 problem objective quadratic and the latter makes it linear. 512 """ 513 )) 514 CONFIG.declare('verbose', ConfigValue( 515 default=False, 516 domain=bool, 517 description="Flag to enable verbose output", 518 doc=""" 519 If True, prints subproblem solutions, as well as potential and added cuts 520 during algorithm. 521 522 If False, only the relaxed BigM objective and minimal information about 523 cuts is logged. 524 """ 525 )) 526 CONFIG.declare('stream_solver', ConfigValue( 527 default=False, 528 domain=bool, 529 description="""If true, sets tee=True for every solve performed over 530 "the course of the algorithm""" 531 )) 532 CONFIG.declare('solver_options', ConfigBlock( 533 implicit=True, 534 description="Dictionary of solver options", 535 doc=""" 536 Dictionary of solver options that will be set for the solver for both the 537 relaxed BigM and separation problem solves. 538 """ 539 )) 540 CONFIG.declare('tighten_relaxation', ConfigValue( 541 default=do_not_tighten, 542 description="Callback which takes the GDP formulation and returns a " 543 "GDP formulation with a tighter hull relaxation", 544 doc=""" 545 Function which accepts the GDP formulation of the problem and returns 546 a GDP formulation which the transformation will then take the hull 547 reformulation of. 548 549 Most typically, this callback would be used to apply basic steps before 550 taking the hull reformulation, but anything which tightens the GDP can 551 be performed here. 552 """ 553 )) 554 CONFIG.declare('create_cuts', ConfigValue( 555 default=create_cuts_normal_vector, 556 description="Callback which generates a list of cuts, given the solved " 557 "relaxed bigM and relaxed hull solutions. If no cuts can be " 558 "generated, returns None", 559 doc=""" 560 Callback to generate cuts to be added to the bigM problem based on 561 solutions to the relaxed bigM problem and the separation problem. 562 563 Arguments 564 --------- 565 transBlock_rBigm: transformation block on relaxed bigM instance 566 transBlock_rHull: transformation blcok on relaxed hull instance 567 var_info: List of tuples (rBigM_var, rHull_var, xstar_param) 568 hull_to_bigm_map: For expression substition, maps id(hull_var) to 569 coresponding bigm var 570 rBigM_linear_constraints: list of linear constraints in relaxed bigM 571 rHull_vars: list of all variables in relaxed hull 572 disaggregated_vars: ComponentSet of disaggregated variables in hull 573 reformulation 574 cut_threshold: Amount x* needs to be infeasible in generated cut in order 575 to consider the cut for addition to the bigM model. 576 zero_tolerance: Tolerance at which a float will be treated as 0 577 578 Returns 579 ------- 580 list of cuts to be added to bigM problem (and relaxed bigM problem), 581 represented as expressions using variables from the bigM model 582 """ 583 )) 584 CONFIG.declare('post_process_cut', ConfigValue( 585 default=back_off_constraint_with_calculated_cut_violation, 586 description="Callback which takes a generated cut and post processes " 587 "it, presumably to back it off in the case of numerical error. Set to " 588 "None if not post-processing is desired.", 589 doc=""" 590 Callback to adjust a cut returned from create_cuts before adding it to 591 the model, presumably to make it more conservative in case of numerical 592 error. 593 594 Arguments 595 --------- 596 cut: the cut to be made more conservative, a Constraint 597 transBlock_rHull: the relaxed hull model's transformation Block. 598 bigm_to_hull_map: Dictionary mapping ids of bigM variables to the 599 corresponding variables on the relaxed hull instance. 600 opt: SolverFactory object for subproblem solves in this procedure 601 stream_solver: Whether or not to set tee=True while solving. 602 TOL: A tolerance 603 604 Returns 605 ------- 606 None, modifies the cut in place 607 """ 608 )) 609 # back off problem tolerance (on top of the solver's (sometimes)) 610 CONFIG.declare('back_off_problem_tolerance', ConfigValue( 611 default=1e-8, 612 domain=NonNegativeFloat, 613 description="Tolerance to pass to the post_process_cut callback.", 614 doc=""" 615 Tolerance passed to the post_process_cut callback. 616 617 Depending on the callback, different values could make sense, but 618 something on the order of the solver's optimality or constraint 619 tolerances is appropriate. 620 """ 621 )) 622 CONFIG.declare('cut_filtering_threshold', ConfigValue( 623 default=0.001, 624 domain=NonNegativeFloat, 625 description="Tolerance used to decide if a cut removes x* from the " 626 "relaxed BigM problem by enough to be added to the bigM problem.", 627 doc=""" 628 Absolute tolerance used to decide whether to keep a cut. We require 629 that, when evaluated at x* (the relaxed BigM optimal solution), the 630 cut be infeasible by at least this tolerance. 631 """ 632 )) 633 CONFIG.declare('zero_tolerance', ConfigValue( 634 default=1e-9, 635 domain=NonNegativeFloat, 636 description="Tolerance at which floats are assumed to be 0 while " 637 "performing Fourier-Motzkin elimination", 638 doc=""" 639 Only relevant when create_cuts=create_cuts_fme, this sets the 640 zero_tolerance option for the Fourier-Motzkin elimination transformation. 641 """ 642 )) 643 CONFIG.declare('do_integer_arithmetic', ConfigValue( 644 default=False, 645 domain=bool, 646 description="Only relevant if using Fourier-Motzkin Elimination (FME) " 647 "and if all problem data is integer, requires FME transformation to " 648 "perform integer arithmetic to eliminate numerical error.", 649 doc=""" 650 Only relevant when create_cuts=create_cuts_fme and if all problem data 651 is integer, this sets the do_integer_arithmetic flag to true for the 652 FME transformation, meaning that the projection to the big-M space 653 can be done with exact precision. 654 """ 655 )) 656 CONFIG.declare('cuts_name', ConfigValue( 657 default=None, 658 domain=str, 659 description="Optional name for the IndexedConstraint containing the " 660 "projected cuts. Must be a unique name with respect to the " 661 "instance.", 662 doc=""" 663 Optional name for the IndexedConstraint containing the projected 664 constraints. If not specified, the cuts will be stored on a 665 private block created by the transformation, so if you want access 666 to them after the transformation, use this argument. 667 668 Must be a string which is a unique component name with respect to the 669 Block on which the transformation is called. 670 """ 671 )) 672 CONFIG.declare('tight_constraint_tolerance', ConfigValue( 673 default=1e-6, # Gurobi constraint tolerance 674 domain=NonNegativeFloat, 675 description="Tolerance at which a constraint is considered tight for " 676 "the Fourier-Motzkin cut generation procedure.", 677 doc=""" 678 For a constraint a^Tx <= b, the Fourier-Motzkin cut generation procedure 679 will consider the constraint tight (and add it to the set of constraints 680 being projected) when a^Tx - b is less than this tolerance. 681 682 It is recommended to set this tolerance to the constraint tolerance of 683 the solver being used. 684 """ 685 )) 686 def __init__(self): 687 super(CuttingPlane_Transformation, self).__init__() 688 689 def _apply_to(self, instance, bigM=None, **kwds): 690 original_log_level = logger.level 691 log_level = logger.getEffectiveLevel() 692 try: 693 assert not NAME_BUFFER 694 self._config = self.CONFIG(kwds.pop('options', {})) 695 self._config.set_value(kwds) 696 697 if self._config.verbose and log_level > logging.INFO: 698 logger.setLevel(logging.INFO) 699 self.verbose = True 700 elif log_level <= logging.INFO: 701 self.verbose = True 702 else: 703 self.verbose = False 704 705 (instance_rBigM, cuts_obj, instance_rHull, var_info, 706 transBlockName) = self._setup_subproblems( instance, bigM, 707 self._config.\ 708 tighten_relaxation) 709 710 self._generate_cuttingplanes( instance_rBigM, cuts_obj, 711 instance_rHull, var_info, 712 transBlockName) 713 714 # restore integrality 715 TransformationFactory('core.relax_integer_vars').apply_to(instance, 716 undo=True) 717 finally: 718 del self._config 719 del self.verbose 720 # clear the global name buffer 721 NAME_BUFFER.clear() 722 # restore logging level 723 logger.setLevel(original_log_level) 724 725 def _setup_subproblems(self, instance, bigM, tighten_relaxation_callback): 726 # create transformation block 727 transBlockName, transBlock = self._add_transformation_block(instance) 728 729 # We store a list of all vars so that we can efficiently 730 # generate maps among the subproblems 731 732 transBlock.all_vars = list(v for v in instance.component_data_objects( 733 Var, 734 descend_into=(Block, Disjunct), 735 sort=SortComponents.deterministic) if not v.is_fixed()) 736 737 # we'll store all the cuts we add together 738 nm = self._config.cuts_name 739 if nm is None: 740 cuts_obj = transBlock.cuts = Constraint(NonNegativeIntegers) 741 else: 742 # check that this really is an available name 743 if instance.component(nm) is not None: 744 raise GDP_Error("cuts_name was specified as '%s', but this is " 745 "already a component on the instance! Please " 746 "specify a unique name." % nm) 747 instance.add_component(nm, Constraint(NonNegativeIntegers)) 748 cuts_obj = instance.component(nm) 749 750 # get bigM and hull relaxations 751 bigMRelaxation = TransformationFactory('gdp.bigm') 752 hullRelaxation = TransformationFactory('gdp.hull') 753 relaxIntegrality = TransformationFactory('core.relax_integer_vars') 754 755 # 756 # Generate the Hull relaxation (used for the separation 757 # problem to generate cutting planes) 758 # 759 tighter_instance = tighten_relaxation_callback(instance) 760 instance_rHull = hullRelaxation.create_using(tighter_instance) 761 relaxIntegrality.apply_to(instance_rHull, 762 transform_deactivated_blocks=True) 763 764 # 765 # Reformulate the instance using the BigM relaxation (this will 766 # be the final instance returned to the user) 767 # 768 bigMRelaxation.apply_to(instance, bigM=bigM) 769 770 # 771 # Generate the continuous relaxation of the BigM transformation. We'll 772 # restore it at the end. 773 # 774 relaxIntegrality.apply_to(instance, transform_deactivated_blocks=True) 775 776 # 777 # Add the xstar parameter for the Hull problem 778 # 779 transBlock_rHull = instance_rHull.component(transBlockName) 780 # 781 # this will hold the solution to rbigm each time we solve it. We 782 # add it to the transformation block so that we don't have to 783 # worry about name conflicts. 784 transBlock_rHull.xstar = Param( range(len(transBlock.all_vars)), 785 mutable=True, default=0, within=Reals) 786 787 # we will add a block that we will deactivate to use to store the 788 # extended space cuts. We never need to solve these, but we need them to 789 # be constructed for the sake of Fourier-Motzkin Elimination 790 extendedSpaceCuts = transBlock_rHull.extendedSpaceCuts = Block() 791 extendedSpaceCuts.deactivate() 792 extendedSpaceCuts.cuts = Constraint(Any) 793 794 # 795 # Generate the mapping between the variables on all the 796 # instances and the xstar parameter. 797 # 798 var_info = [ 799 (v, # this is the bigM variable 800 transBlock_rHull.all_vars[i], 801 transBlock_rHull.xstar[i]) 802 for i,v in enumerate(transBlock.all_vars)] 803 804 # NOTE: we wait to add the separation objective to the rHull problem 805 # because it is best to do it in the first iteration, so that we can 806 # skip stale variables. 807 808 return (instance, cuts_obj, instance_rHull, var_info, transBlockName) 809 810 # this is the map that I need to translate my projected cuts and add 811 # them to bigM 812 def _create_hull_to_bigm_substitution_map(self, var_info): 813 return dict((id(var_info[i][1]), var_info[i][0]) for i in 814 range(len(var_info))) 815 816 # this map is needed to solve the back-off problem for post-processing 817 def _create_bigm_to_hull_substition_map(self, var_info): 818 return dict((id(var_info[i][0]), var_info[i][1]) for i in 819 range(len(var_info))) 820 821 def _get_disaggregated_vars(self, hull): 822 disaggregatedVars = ComponentSet() 823 hull_xform = TransformationFactory('gdp.hull') 824 for disjunction in hull.component_data_objects( Disjunction, 825 descend_into=(Disjunct, 826 Block)): 827 for disjunct in disjunction.disjuncts: 828 if disjunct.transformation_block is not None: 829 transBlock = disjunct.transformation_block() 830 for v in transBlock.disaggregatedVars.\ 831 component_data_objects(Var): 832 disaggregatedVars.add(v) 833 834 return disaggregatedVars 835 836 def _get_rBigM_obj_and_constraints(self, instance_rBigM): 837 # We try to grab the first active objective. If there is more 838 # than one, the writer will yell when we try to solve below. If 839 # there are 0, we will yell here. 840 rBigM_obj = next(instance_rBigM.component_data_objects( 841 Objective, active=True), None) 842 if rBigM_obj is None: 843 raise GDP_Error("Cannot apply cutting planes transformation " 844 "without an active objective in the model!") 845 846 # 847 # Collect all of the linear constraints that are in the rBigM 848 # instance. We will need these so that we can compare what we get from 849 # FME to them and make sure we aren't adding redundant constraints to 850 # the model. For convenience, we will make sure they are all in the form 851 # lb <= expr (so we will break equality constraints) 852 # 853 fme = TransformationFactory('contrib.fourier_motzkin_elimination') 854 rBigM_linear_constraints = [] 855 for cons in instance_rBigM.component_data_objects( 856 Constraint, 857 descend_into=Block, 858 sort=SortComponents.deterministic, 859 active=True): 860 body = cons.body 861 if body.polynomial_degree() != 1: 862 # We will never get a nonlinear constraint out of FME, so we 863 # don't risk it being identical to this one. 864 continue 865 866 # TODO: Guess this shouldn't have been private... 867 rBigM_linear_constraints.extend(fme._process_constraint(cons)) 868 869 # [ESJ Aug 13 2020] NOTE: We actually don't need to worry about variable 870 # bounds here because the FME transformation will take care of them 871 # (i.e. convert those of the disaggregated variables to constraints for 872 # the purposes of the projection.) 873 874 return rBigM_obj, rBigM_linear_constraints 875 876 def _generate_cuttingplanes( self, instance_rBigM, cuts_obj, instance_rHull, 877 var_info, transBlockName): 878 879 opt = SolverFactory(self._config.solver) 880 stream_solver = self._config.stream_solver 881 opt.options = dict(self._config.solver_options) 882 883 improving = True 884 prev_obj = None 885 epsilon = self._config.minimum_improvement_threshold 886 cuts = None 887 888 transBlock_rHull = instance_rHull.component(transBlockName) 889 890 rBigM_obj, rBigM_linear_constraints = self.\ 891 _get_rBigM_obj_and_constraints( 892 instance_rBigM) 893 894 # Get list of all variables in the rHull model which we will use when 895 # calculating the composite normal vector. 896 rHull_vars = [i for i in instance_rHull.component_data_objects( 897 Var, 898 descend_into=Block, 899 sort=SortComponents.deterministic)] 900 901 # collect a list of disaggregated variables. 902 disaggregated_vars = self._get_disaggregated_vars( instance_rHull) 903 904 hull_to_bigm_map = self._create_hull_to_bigm_substitution_map(var_info) 905 bigm_to_hull_map = self._create_bigm_to_hull_substition_map(var_info) 906 xhat = ComponentMap() 907 908 while (improving): 909 # solve rBigM, solution is xstar 910 results = opt.solve(instance_rBigM, tee=stream_solver, 911 load_solutions=False) 912 if verify_successful_solve(results) is not NORMAL: 913 logger.warning("Relaxed BigM subproblem " 914 "did not solve normally. Stopping cutting " 915 "plane generation.\n\n%s" % (results,)) 916 return 917 instance_rBigM.solutions.load_from(results) 918 919 rBigM_objVal = value(rBigM_obj) 920 logger.warning("rBigM objective = %s" % (rBigM_objVal,)) 921 922 # 923 # Add the separation objective to the hull subproblem if it's not 924 # there already (so in the first iteration). We're waiting until now 925 # to avoid it including variables that came back stale from the 926 # rbigm solve. 927 # 928 if transBlock_rHull.component("separation_objective") is None: 929 self._add_separation_objective(var_info, transBlock_rHull) 930 931 # copy over xstar 932 logger.info("x* is:") 933 for x_rbigm, x_hull, x_star in var_info: 934 if not x_rbigm.stale: 935 x_star.value = x_rbigm.value 936 # initialize the X values 937 x_hull.value = x_rbigm.value 938 if self.verbose: 939 logger.info("\t%s = %s" % 940 (x_rbigm.getname(fully_qualified=True, 941 name_buffer=NAME_BUFFER), 942 x_rbigm.value)) 943 944 # compare objectives: check absolute difference close to 0, relative 945 # difference further from 0. 946 if prev_obj is None: 947 improving = True 948 else: 949 obj_diff = prev_obj - rBigM_objVal 950 improving = ( abs(obj_diff) > epsilon if abs(rBigM_objVal) < 1 951 else abs(obj_diff/prev_obj) > epsilon ) 952 953 # solve separation problem to get xhat. 954 results = opt.solve(instance_rHull, tee=stream_solver, 955 load_solutions=False) 956 if verify_successful_solve(results) is not NORMAL: 957 logger.warning("Hull separation subproblem " 958 "did not solve normally. Stopping cutting " 959 "plane generation.\n\n%s" % (results,)) 960 return 961 instance_rHull.solutions.load_from(results) 962 logger.warning("separation problem objective value: %s" % 963 value(transBlock_rHull.separation_objective)) 964 965 # save xhat to initialize rBigM with in the next iteration 966 if self.verbose: 967 logger.info("xhat is: ") 968 for x_rbigm, x_hull, x_star in var_info: 969 xhat[x_rbigm] = value(x_hull) 970 if self.verbose: 971 logger.info("\t%s = %s" % 972 (x_hull.getname(fully_qualified=True, 973 name_buffer=NAME_BUFFER), 974 x_hull.value)) 975 976 # [JDS 19 Dec 18] Note: we check that the separation objective was 977 # significantly nonzero. If it is too close to zero, either the 978 # rBigM solution was in the convex hull, or the separation vector is 979 # so close to zero that the resulting cut is likely to have 980 # numerical issues. 981 if value(transBlock_rHull.separation_objective) < \ 982 self._config.separation_objective_threshold: 983 logger.warning("Separation problem objective below threshold of" 984 " %s: Stopping cut generation." % 985 self._config.separation_objective_threshold) 986 break 987 988 cuts = self._config.create_cuts(transBlock_rHull, var_info, 989 hull_to_bigm_map, 990 rBigM_linear_constraints, 991 rHull_vars, disaggregated_vars, 992 self._config.norm, 993 self._config.cut_filtering_threshold, 994 self._config.zero_tolerance, 995 self._config.do_integer_arithmetic, 996 self._config.\ 997 tight_constraint_tolerance) 998 999 # We are done if the cut generator couldn't return a valid cut 1000 if cuts is None: 1001 logger.warning("Did not generate a valid cut, stopping cut " 1002 "generation.") 1003 break 1004 if not improving: 1005 logger.warning("Difference in relaxed BigM problem objective " 1006 "values from past two iterations is below " 1007 "threshold of %s: Stopping cut generation." % 1008 epsilon) 1009 break 1010 1011 for cut in cuts: 1012 # we add the cut to the model and then post-process it in place. 1013 cut_number = len(cuts_obj) 1014 logger.warning("Adding cut %s to BigM model." % (cut_number,)) 1015 cuts_obj.add(cut_number, cut) 1016 if self._config.post_process_cut is not None: 1017 self._config.post_process_cut( 1018 cuts_obj[cut_number], transBlock_rHull, 1019 bigm_to_hull_map, opt, stream_solver, 1020 self._config.back_off_problem_tolerance) 1021 1022 if cut_number + 1 == self._config.max_number_of_cuts: 1023 logger.warning("Reached maximum number of cuts.") 1024 break 1025 1026 prev_obj = rBigM_objVal 1027 1028 # Initialize rbigm with xhat (for the next iteration) 1029 for x_rbigm, x_hull, x_star in var_info: 1030 x_rbigm.value = xhat[x_rbigm] 1031 1032 def _add_transformation_block(self, instance): 1033 # creates transformation block with a unique name based on name, adds it 1034 # to instance, and returns it. 1035 transBlockName = unique_component_name( 1036 instance, 1037 '_pyomo_gdp_cuttingplane_transformation') 1038 transBlock = Block() 1039 instance.add_component(transBlockName, transBlock) 1040 return transBlockName, transBlock 1041 1042 1043 def _add_separation_objective(self, var_info, transBlock_rHull): 1044 # creates the separation objective. That is just adding an objective for 1045 # Euclidean norm, it means adding an auxilary variable to linearize the 1046 # L-infinity norm. We do this assuming that rBigM has been solved, and 1047 # if any variables come back stale, we leave them out of the separation 1048 # problem, as they aren't doing anything and they could cause numerical 1049 # issues later. 1050 1051 # Deactivate any/all other objectives 1052 for o in transBlock_rHull.model().component_data_objects(Objective): 1053 o.deactivate() 1054 norm = self._config.norm 1055 to_delete = [] 1056 1057 if norm == 2: 1058 obj_expr = 0 1059 for i, (x_rbigm, x_hull, x_star) in enumerate(var_info): 1060 if not x_rbigm.stale: 1061 obj_expr += (x_hull - x_star)**2 1062 else: 1063 if self.verbose: 1064 logger.info("The variable %s will not be included in " 1065 "the separation problem: It was stale in " 1066 "the rBigM solve." % x_rbigm.getname( 1067 fully_qualified=True, 1068 name_buffer=NAME_BUFFER)) 1069 to_delete.append(i) 1070 elif norm == float('inf'): 1071 u = transBlock_rHull.u = Var(domain=NonNegativeReals) 1072 inf_cons = transBlock_rHull.inf_norm_linearization = Constraint( 1073 NonNegativeIntegers) 1074 i = 0 1075 for j, (x_rbigm, x_hull, x_star) in enumerate(var_info): 1076 if not x_rbigm.stale: 1077 # NOTE: these are written as >= constraints so that we know 1078 # the duals will come back nonnegative. 1079 inf_cons[i] = u - x_hull >= - x_star 1080 inf_cons[i+1] = u + x_hull >= x_star 1081 i += 2 1082 else: 1083 if self.verbose: 1084 logger.info("The variable %s will not be included in " 1085 "the separation problem: It was stale in " 1086 "the rBigM solve." % x_rbigm.getname( 1087 fully_qualified=True, 1088 name_buffer=NAME_BUFFER)) 1089 to_delete.append(j) 1090 # we'll need the duals of these to get the subgradient 1091 self._add_dual_suffix(transBlock_rHull.model()) 1092 obj_expr = u 1093 1094 # delete the unneeded x_stars so that we don't add cuts involving 1095 # useless variables later. 1096 for i in sorted(to_delete, reverse=True): 1097 del var_info[i] 1098 1099 # add separation objective to transformation block 1100 transBlock_rHull.separation_objective = Objective(expr=obj_expr) 1101 1102 def _add_dual_suffix(self, rHull): 1103 # rHull is our model and we aren't giving it back (unless in the future 1104 # we we add a callback to do basic steps to it...), so we just check if 1105 # dual is there. If it's a Suffix, we'll borrow it. If it's something 1106 # else we'll rename it and add the Suffix. 1107 dual = rHull.component("dual") 1108 if dual is None: 1109 rHull.dual = Suffix(direction=Suffix.IMPORT) 1110 else: 1111 if dual.ctype is Suffix: 1112 return 1113 rHull.del_component(dual) 1114 rHull.dual = Suffix(direction=Suffix.IMPORT) 1115 rHull.add_component(unique_component_name(rHull, "dual"), dual) 1116