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"""Utility functions and classes for the GDPopt solver.""" 11 12from __future__ import division 13 14import logging 15from contextlib import contextmanager 16from math import fabs 17 18from pyomo.common import deprecated, timing 19from pyomo.common.collections import ComponentSet, Bunch 20from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr 21from pyomo.contrib.gdpopt.data_class import GDPoptSolveData 22from pyomo.contrib.mcpp.pyomo_mcpp import mcpp_available, McCormick 23from pyomo.core import (Block, Constraint, 24 Objective, Reals, Var, minimize, value, ConstraintList) 25from pyomo.core.expr.current import identify_variables 26from pyomo.gdp import Disjunct, Disjunction 27from pyomo.opt import SolverFactory, SolverResults 28from pyomo.opt.results import ProblemSense 29from pyomo.util.model_size import build_model_size_report 30 31 32class _DoNothing(object): 33 """Do nothing, literally. 34 35 This class is used in situations of "do something if attribute exists." 36 """ 37 38 def __init__(self, *args, **kwargs): 39 pass 40 41 def __call__(self, *args, **kwargs): 42 pass 43 44 def __getattr__(self, attr): 45 def _do_nothing(*args, **kwargs): 46 pass 47 48 return _do_nothing 49 50 51class SuppressInfeasibleWarning(object): 52 """Suppress the infeasible model warning message from solve(). 53 54 The "WARNING: Loading a SolverResults object with a warning status" warning 55 message from calling solve() is often unwanted, but there is no clear way 56 to suppress it. 57 58 This is modeled on LoggingIntercept from pyomo.common.log, 59 but different in function. 60 61 """ 62 63 class InfeasibleWarningFilter(logging.Filter): 64 def filter(self, record): 65 return not record.getMessage().startswith( 66 "Loading a SolverResults object with a warning status into model=") 67 68 warning_filter = InfeasibleWarningFilter() 69 70 def __enter__(self): 71 logger = logging.getLogger('pyomo.core') 72 logger.addFilter(self.warning_filter) 73 74 def __exit__(self, exception_type, exception_value, traceback): 75 logger = logging.getLogger('pyomo.core') 76 logger.removeFilter(self.warning_filter) 77 78 79def presolve_lp_nlp(solve_data, config): 80 """If the model is an LP or NLP, solve it directly. 81 82 """ 83 m = solve_data.working_model 84 GDPopt = m.GDPopt_utils 85 86 # Handle LP/NLP being passed to the solver 87 prob = solve_data.results.problem 88 if (prob.number_of_binary_variables == 0 and 89 prob.number_of_integer_variables == 0 and 90 prob.number_of_disjunctions == 0): 91 config.logger.info('Problem has no discrete decisions.') 92 obj = next(m.component_data_objects(Objective, active=True)) 93 if (any(c.body.polynomial_degree() not in (1, 0) for c in GDPopt.constraint_list) or 94 obj.expr.polynomial_degree() not in (1, 0)): 95 config.logger.info( 96 "Your model is an NLP (nonlinear program). " 97 "Using NLP solver %s to solve." % config.nlp_solver) 98 results = SolverFactory(config.nlp_solver).solve( 99 solve_data.original_model, **config.nlp_solver_args) 100 return True, results 101 else: 102 config.logger.info( 103 "Your model is an LP (linear program). " 104 "Using LP solver %s to solve." % config.mip_solver) 105 results = SolverFactory(config.mip_solver).solve( 106 solve_data.original_model, **config.mip_solver_args) 107 return True, results 108 109 # TODO if any continuous variables are multipled with binary ones, need 110 # to do some kind of transformation (Glover?) or throw an error message 111 return False, None 112 113 114def process_objective(solve_data, config, move_linear_objective=False, use_mcpp=True, updata_var_con_list=True): 115 """Process model objective function. 116 117 Check that the model has only 1 valid objective. 118 If the objective is nonlinear, move it into the constraints. 119 If no objective function exists, emit a warning and create a dummy objective. 120 121 Parameters 122 ---------- 123 solve_data (GDPoptSolveData): solver environment data class 124 config (ConfigBlock): solver configuration options 125 move_linear_objective (bool): if True, move even linear 126 objective functions to the constraints 127 128 """ 129 m = solve_data.working_model 130 util_blk = getattr(m, solve_data.util_block_name) 131 # Handle missing or multiple objectives 132 active_objectives = list(m.component_data_objects( 133 ctype=Objective, active=True, descend_into=True)) 134 solve_data.results.problem.number_of_objectives = len(active_objectives) 135 if len(active_objectives) == 0: 136 config.logger.warning( 137 'Model has no active objectives. Adding dummy objective.') 138 util_blk.dummy_objective = Objective(expr=1) 139 main_obj = util_blk.dummy_objective 140 elif len(active_objectives) > 1: 141 raise ValueError('Model has multiple active objectives.') 142 else: 143 main_obj = active_objectives[0] 144 solve_data.results.problem.sense = ProblemSense.minimize if main_obj.sense == 1 else ProblemSense.maximize 145 solve_data.objective_sense = main_obj.sense 146 147 # Move the objective to the constraints if it is nonlinear 148 if main_obj.expr.polynomial_degree() not in (1, 0) \ 149 or move_linear_objective: 150 if move_linear_objective: 151 config.logger.info("Moving objective to constraint set.") 152 else: 153 config.logger.info( 154 "Objective is nonlinear. Moving it to constraint set.") 155 156 util_blk.objective_value = Var(domain=Reals, initialize=0) 157 if mcpp_available() and use_mcpp: 158 mc_obj = McCormick(main_obj.expr) 159 util_blk.objective_value.setub(mc_obj.upper()) 160 util_blk.objective_value.setlb(mc_obj.lower()) 161 else: 162 # Use Pyomo's contrib.fbbt package 163 lb, ub = compute_bounds_on_expr(main_obj.expr) 164 if solve_data.results.problem.sense == ProblemSense.minimize: 165 util_blk.objective_value.setlb(lb) 166 else: 167 util_blk.objective_value.setub(ub) 168 169 if main_obj.sense == minimize: 170 util_blk.objective_constr = Constraint( 171 expr=util_blk.objective_value >= main_obj.expr) 172 else: 173 util_blk.objective_constr = Constraint( 174 expr=util_blk.objective_value <= main_obj.expr) 175 # Deactivate the original objective and add this new one. 176 main_obj.deactivate() 177 util_blk.objective = Objective( 178 expr=util_blk.objective_value, sense=main_obj.sense) 179 # Add the new variable and constraint to the working lists 180 if main_obj.expr.polynomial_degree() not in (1, 0) or (move_linear_objective and updata_var_con_list): 181 util_blk.variable_list.append(util_blk.objective_value) 182 util_blk.continuous_variable_list.append(util_blk.objective_value) 183 util_blk.constraint_list.append(util_blk.objective_constr) 184 util_blk.objective_list.append(util_blk.objective) 185 if util_blk.objective_constr.body.polynomial_degree() in (0, 1): 186 util_blk.linear_constraint_list.append(util_blk.objective_constr) 187 else: 188 util_blk.nonlinear_constraint_list.append( 189 util_blk.objective_constr) 190 191 192def a_logger(str_or_logger): 193 """Returns a logger when passed either a logger name or logger object.""" 194 if isinstance(str_or_logger, logging.Logger): 195 return str_or_logger 196 else: 197 return logging.getLogger(str_or_logger) 198 199 200def copy_var_list_values(from_list, to_list, config, 201 skip_stale=False, skip_fixed=True, 202 ignore_integrality=False): 203 """Copy variable values from one list to another. 204 205 Rounds to Binary/Integer if neccessary 206 Sets to zero for NonNegativeReals if neccessary 207 """ 208 for v_from, v_to in zip(from_list, to_list): 209 if skip_stale and v_from.stale: 210 continue # Skip stale variable values. 211 if skip_fixed and v_to.is_fixed(): 212 continue # Skip fixed variables. 213 try: 214 v_to.set_value(value(v_from, exception=False)) 215 if skip_stale: 216 v_to.stale = False 217 except ValueError as err: 218 err_msg = getattr(err, 'message', str(err)) 219 var_val = value(v_from) 220 rounded_val = int(round(var_val)) 221 # Check to see if this is just a tolerance issue 222 if ignore_integrality \ 223 and v_to.is_integer(): # not v_to.is_continuous() 224 v_to.value = value(v_from, exception=False) 225 elif v_to.is_integer() and (fabs(var_val - rounded_val) <= config.integer_tolerance): # not v_to.is_continuous() 226 v_to.set_value(rounded_val) 227 elif 'is not in domain NonNegativeReals' in err_msg and ( 228 fabs(var_val) <= config.zero_tolerance): 229 v_to.set_value(0) 230 else: 231 raise 232 233 234def is_feasible(model, config): 235 """Checks to see if the algebraic model is feasible in its current state. 236 237 Checks variable bounds and active constraints. Not for use with 238 untransformed GDP models. 239 240 """ 241 disj = next(model.component_data_objects( 242 ctype=Disjunct, active=True), None) 243 if disj is not None: 244 raise NotImplementedError( 245 "Found active disjunct %s. " 246 "This function is not intended to check " 247 "feasibility of disjunctive models, " 248 "only transformed subproblems." % disj.name) 249 250 config.logger.debug('Checking if model is feasible.') 251 for constr in model.component_data_objects( 252 ctype=Constraint, active=True, descend_into=True): 253 # Check constraint lower bound 254 if (constr.lower is not None and ( 255 value(constr.lower) - value(constr.body) 256 >= config.constraint_tolerance 257 )): 258 config.logger.info('%s: body %s < LB %s' % ( 259 constr.name, value(constr.body), value(constr.lower))) 260 return False 261 # check constraint upper bound 262 if (constr.upper is not None and ( 263 value(constr.body) - value(constr.upper) 264 >= config.constraint_tolerance 265 )): 266 config.logger.info('%s: body %s > UB %s' % ( 267 constr.name, value(constr.body), value(constr.upper))) 268 return False 269 for var in model.component_data_objects(ctype=Var, descend_into=True): 270 # Check variable lower bound 271 if (var.has_lb() and 272 value(var.lb) - value(var) >= config.variable_tolerance): 273 config.logger.info('%s: value %s < LB %s' % ( 274 var.name, value(var), value(var.lb))) 275 return False 276 # Check variable upper bound 277 if (var.has_ub() and 278 value(var) - value(var.ub) >= config.variable_tolerance): 279 config.logger.info('%s: value %s > UB %s' % ( 280 var.name, value(var), value(var.ub))) 281 return False 282 config.logger.info('Model is feasible.') 283 return True 284 285 286def build_ordered_component_lists(model, solve_data): 287 """Define lists used for future data transfer. 288 289 Also attaches ordered lists of the variables, constraints, disjuncts, and 290 disjunctions to the model so that they can be used for mapping back and 291 forth. 292 293 """ 294 util_blk = getattr(model, solve_data.util_block_name) 295 var_set = ComponentSet() 296 setattr( 297 util_blk, 'constraint_list', list( 298 model.component_data_objects( 299 ctype=Constraint, active=True, 300 descend_into=(Block, Disjunct)))) 301 # print(util_blk.constraint_list) 302 setattr( 303 util_blk, 'linear_constraint_list', list(c for c in model.component_data_objects( 304 ctype=Constraint, active=True, descend_into=(Block, Disjunct)) 305 if c.body.polynomial_degree() in (0, 1))) 306 # print(util_blk.linear_constraint_list) 307 setattr( 308 util_blk, 'nonlinear_constraint_list', list(c for c in model.component_data_objects( 309 ctype=Constraint, active=True, descend_into=(Block, Disjunct)) 310 if c.body.polynomial_degree() not in (0, 1))) 311 # print(util_blk.nonlinear_constraint_list) 312 setattr( 313 util_blk, 'disjunct_list', list( 314 model.component_data_objects( 315 ctype=Disjunct, active=True, 316 descend_into=(Block, Disjunct)))) 317 setattr( 318 util_blk, 'disjunction_list', list( 319 model.component_data_objects( 320 ctype=Disjunction, active=True, 321 descend_into=(Disjunct, Block)))) 322 setattr( 323 util_blk, 'objective_list', list( 324 model.component_data_objects( 325 ctype=Objective, active=True, 326 descend_into=(Block)))) 327 328 # Identify the non-fixed variables in (potentially) active constraints and 329 # objective functions 330 for constr in getattr(util_blk, 'constraint_list'): 331 for v in identify_variables(constr.body, include_fixed=False): 332 var_set.add(v) 333 for obj in model.component_data_objects(ctype=Objective, active=True): 334 for v in identify_variables(obj.expr, include_fixed=False): 335 var_set.add(v) 336 # Disjunct indicator variables might not appear in active constraints. In 337 # fact, if we consider them Logical variables, they should not appear in 338 # active algebraic constraints. For now, they need to be added to the 339 # variable set. 340 for disj in getattr(util_blk, 'disjunct_list'): 341 var_set.add(disj.binary_indicator_var) 342 343 # We use component_data_objects rather than list(var_set) in order to 344 # preserve a deterministic ordering. 345 var_list = list( 346 v for v in model.component_data_objects( 347 ctype=Var, descend_into=(Block, Disjunct)) 348 if v in var_set) 349 setattr(util_blk, 'variable_list', var_list) 350 discrete_variable_list = list( 351 v for v in model.component_data_objects( 352 ctype=Var, descend_into=(Block, Disjunct)) 353 if v in var_set and v.is_integer()) 354 setattr(util_blk, 'discrete_variable_list', discrete_variable_list) 355 continuous_variable_list = list( 356 v for v in model.component_data_objects( 357 ctype=Var, descend_into=(Block, Disjunct)) 358 if v in var_set and v.is_continuous()) 359 setattr(util_blk, 'continuous_variable_list', continuous_variable_list) 360 361 362def setup_results_object(solve_data, config): 363 """Record problem statistics for original model.""" 364 # Create the solver results object 365 res = solve_data.results 366 prob = res.problem 367 res.problem.name = solve_data.original_model.name 368 res.problem.number_of_nonzeros = None # TODO 369 # TODO work on termination condition and message 370 res.solver.termination_condition = None 371 res.solver.message = None 372 res.solver.user_time = None 373 res.solver.system_time = None 374 res.solver.wallclock_time = None 375 res.solver.termination_message = None 376 377 num_of = build_model_size_report(solve_data.original_model) 378 379 # Get count of constraints and variables 380 prob.number_of_constraints = num_of.activated.constraints 381 prob.number_of_disjunctions = num_of.activated.disjunctions 382 prob.number_of_variables = num_of.activated.variables 383 prob.number_of_binary_variables = num_of.activated.binary_variables 384 prob.number_of_continuous_variables = num_of.activated.continuous_variables 385 prob.number_of_integer_variables = num_of.activated.integer_variables 386 387 config.logger.info( 388 "Original model has %s constraints (%s nonlinear) " 389 "and %s disjunctions, " 390 "with %s variables, of which %s are binary, %s are integer, " 391 "and %s are continuous." % 392 (num_of.activated.constraints, 393 num_of.activated.nonlinear_constraints, 394 num_of.activated.disjunctions, 395 num_of.activated.variables, 396 num_of.activated.binary_variables, 397 num_of.activated.integer_variables, 398 num_of.activated.continuous_variables)) 399 400 401# def validate_disjunctions(model, config): 402# """Validate that the active disjunctions on the model are satisfied 403# by the current disjunct indicator_var values.""" 404# active_disjunctions = model.component_data_objects( 405# ctype=Disjunction, active=True, descend_into=(Block, Disjunct)) 406# for disjtn in active_disjunctions: 407# sum_disj_vals = sum(disj.indicator_var.value 408# for disj in disjtn.disjuncts) 409# if disjtn.xor and fabs(sum_disj_vals - 1) > config.integer_tolerance: 410# raise ValueError( 411# "Expected disjunct values to add up to 1 " 412# "for XOR disjunction %s. " 413# "Instead, values add up to %s." % (disjtn.name, sum_disj_vals)) 414# elif sum_disj_vals + config.integer_tolerance < 1: 415# raise ValueError( 416# "Expected disjunct values to add up to at least 1 for " 417# "OR disjunction %s. " 418# "Instead, values add up to %s." % (disjtn.name, sum_disj_vals)) 419 420 421def constraints_in_True_disjuncts(model, config): 422 """Yield constraints in disjuncts where the indicator value is set or fixed to True.""" 423 for constr in model.component_data_objects(Constraint): 424 yield constr 425 observed_disjuncts = ComponentSet() 426 for disjctn in model.component_data_objects(Disjunction): 427 # get all the disjuncts in the disjunction. Check which ones are True. 428 for disj in disjctn.disjuncts: 429 if disj in observed_disjuncts: 430 continue 431 observed_disjuncts.add(disj) 432 if fabs(disj.binary_indicator_var.value - 1) \ 433 <= config.integer_tolerance: 434 for constr in disj.component_data_objects(Constraint): 435 yield constr 436 437 438@contextmanager 439def time_code(timing_data_obj, code_block_name, is_main_timer=False): 440 """Starts timer at entry, stores elapsed time at exit 441 442 If `is_main_timer=True`, the start time is stored in the timing_data_obj, 443 allowing calculation of total elapsed time 'on the fly' (e.g. to enforce 444 a time limit) using `get_main_elapsed_time(timing_data_obj)`. 445 """ 446 start_time = timing.default_timer() 447 if is_main_timer: 448 timing_data_obj.main_timer_start_time = start_time 449 yield 450 elapsed_time = timing.default_timer() - start_time 451 prev_time = timing_data_obj.get(code_block_name, 0) 452 timing_data_obj[code_block_name] = prev_time + elapsed_time 453 454 455def get_main_elapsed_time(timing_data_obj): 456 """Returns the time since entering the main `time_code` context""" 457 current_time = timing.default_timer() 458 try: 459 return current_time - timing_data_obj.main_timer_start_time 460 except AttributeError as e: 461 if 'main_timer_start_time' in str(e): 462 raise e from AttributeError( 463 "You need to be in a 'time_code' context to use `get_main_elapsed_time()`." 464 ) 465 466 467@deprecated( 468 "'restore_logger_level()' has been deprecated in favor of the more " 469 "specific 'lower_logger_level_to()' function.", 470 version='5.6.9') 471@contextmanager 472def restore_logger_level(logger): 473 old_logger_level = logger.level 474 yield 475 logger.setLevel(old_logger_level) 476 477 478@contextmanager 479def lower_logger_level_to(logger, level=None): 480 """Increases logger verbosity by lowering reporting level.""" 481 if level is not None and logger.getEffectiveLevel() > level: 482 # If logger level is higher (less verbose), decrease it 483 old_logger_level = logger.level 484 logger.setLevel(level) 485 yield 486 logger.setLevel(old_logger_level) 487 else: 488 yield # Otherwise, leave the logger alone 489 490 491@contextmanager 492def create_utility_block(model, name, solve_data): 493 created_util_block = False 494 # Create a model block on which to store GDPopt-specific utility 495 # modeling objects. 496 if hasattr(model, name): 497 raise RuntimeError( 498 "GDPopt needs to create a Block named %s " 499 "on the model object, but an attribute with that name " 500 "already exists." % name) 501 else: 502 created_util_block = True 503 setattr(model, name, Block( 504 doc="Container for GDPopt solver utility modeling objects")) 505 solve_data.util_block_name = name 506 507 # Save ordered lists of main modeling components, so that data can 508 # be easily transferred between future model clones. 509 build_ordered_component_lists(model, solve_data) 510 yield 511 if created_util_block: 512 model.del_component(name) 513 514 515@contextmanager 516def setup_solver_environment(model, config): 517 solve_data = GDPoptSolveData() # data object for storing solver state 518 solve_data.config = config 519 solve_data.results = SolverResults() 520 solve_data.timing = Bunch() 521 min_logging_level = logging.INFO if config.tee else None 522 with time_code(solve_data.timing, 'total', is_main_timer=True), \ 523 lower_logger_level_to(config.logger, min_logging_level), \ 524 create_utility_block(model, 'GDPopt_utils', solve_data): 525 526 # Create a working copy of the original model 527 solve_data.original_model = model 528 solve_data.working_model = model.clone() 529 setup_results_object(solve_data, config) 530 solve_data.active_strategy = config.strategy 531 util_block = solve_data.working_model.GDPopt_utils 532 533 # Save model initial values. 534 # These can be used later to initialize NLP subproblems. 535 solve_data.initial_var_values = list( 536 v.value for v in util_block.variable_list) 537 solve_data.best_solution_found = None 538 539 # Integer cuts exclude particular discrete decisions 540 util_block.integer_cuts = ConstraintList(doc='integer cuts') 541 542 # Set up iteration counters 543 solve_data.master_iteration = 0 544 solve_data.mip_iteration = 0 545 solve_data.nlp_iteration = 0 546 547 # set up bounds 548 solve_data.LB = float('-inf') 549 solve_data.UB = float('inf') 550 solve_data.iteration_log = {} 551 552 # Flag indicating whether the solution improved in the past 553 # iteration or not 554 solve_data.feasible_solution_improved = False 555 556 yield solve_data # yield setup solver environment 557 558 if (solve_data.best_solution_found is not None 559 and solve_data.best_solution_found is not solve_data.original_model): 560 # Update values on the original model 561 copy_var_list_values( 562 from_list=solve_data.best_solution_found.GDPopt_utils.variable_list, 563 to_list=solve_data.original_model.GDPopt_utils.variable_list, 564 config=config) 565 566 # Finalize results object 567 solve_data.results.problem.lower_bound = solve_data.LB 568 solve_data.results.problem.upper_bound = solve_data.UB 569 solve_data.results.solver.iterations = solve_data.master_iteration 570 solve_data.results.solver.timing = solve_data.timing 571 solve_data.results.solver.user_time = solve_data.timing.total 572 solve_data.results.solver.wallclock_time = solve_data.timing.total 573 574 575def indent(text, prefix): 576 """This should be replaced with textwrap.indent when we stop supporting python 2.7.""" 577 return ''.join(prefix + line for line in text.splitlines(True)) 578