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"""main problem functions.""" 12from __future__ import division 13import logging 14from pyomo.core import Constraint, Expression, Objective, minimize, value 15from pyomo.opt import TerminationCondition as tc 16from pyomo.opt import SolutionStatus, SolverFactory 17from pyomo.contrib.gdpopt.util import copy_var_list_values, SuppressInfeasibleWarning, _DoNothing, get_main_elapsed_time, time_code 18from pyomo.contrib.gdpopt.mip_solve import distinguish_mip_infeasible_or_unbounded 19from pyomo.solvers.plugins.solvers.persistent_solver import PersistentSolver 20from pyomo.common.dependencies import attempt_import 21from pyomo.contrib.mindtpy.util import generate_norm1_objective_function, generate_norm2sq_objective_function, generate_norm_inf_objective_function, generate_lag_objective_function, set_solver_options, GurobiPersistent4MindtPy 22 23 24logger = logging.getLogger('pyomo.contrib.mindtpy') 25 26 27single_tree, single_tree_available = attempt_import( 28 'pyomo.contrib.mindtpy.single_tree') 29tabu_list, tabu_list_available = attempt_import( 30 'pyomo.contrib.mindtpy.tabu_list') 31 32 33def solve_main(solve_data, config, fp=False, regularization_problem=False): 34 """ 35 This function solves the MIP main problem 36 37 Parameters 38 ---------- 39 solve_data: MindtPy Data Container 40 data container that holds solve-instance data 41 config: ConfigBlock 42 contains the specific configurations for the algorithm 43 44 Returns 45 ------- 46 solve_data.mip: Pyomo model 47 the MIP stored in solve_data 48 main_mip_results: Pyomo results object 49 result from solving the main MIP 50 fp: Bool 51 generate the feasibility pump regularization main problem 52 regularization_problem: Bool 53 generate the ROA regularization main problem 54 """ 55 if fp: 56 config.logger.info('FP-MIP %s: Solve main problem.' % 57 (solve_data.fp_iter,)) 58 elif regularization_problem: 59 config.logger.info('Regularization-MIP %s: Solve main regularization problem.' % 60 (solve_data.mip_iter,)) 61 else: 62 solve_data.mip_iter += 1 63 config.logger.info('MIP %s: Solve main problem.' % 64 (solve_data.mip_iter,)) 65 66 # setup main problem 67 setup_main(solve_data, config, fp, regularization_problem) 68 mainopt = setup_mip_solver(solve_data, config, regularization_problem) 69 70 mip_args = dict(config.mip_solver_args) 71 if config.mip_solver in {'cplex', 'cplex_persistent', 'gurobi', 'gurobi_persistent'}: 72 mip_args['warmstart'] = True 73 set_solver_options(mainopt, solve_data, config, 74 solver_type='mip', regularization=regularization_problem) 75 try: 76 with time_code(solve_data.timing, 'regularization main' if regularization_problem else ('fp main' if fp else 'main')): 77 main_mip_results = mainopt.solve(solve_data.mip, 78 tee=config.mip_solver_tee, **mip_args) 79 except (ValueError, AttributeError): 80 if config.single_tree: 81 config.logger.warning('Single tree terminate.') 82 if get_main_elapsed_time(solve_data.timing) >= config.time_limit - 2: 83 config.logger.warning('due to the timelimit.') 84 solve_data.results.solver.termination_condition = tc.maxTimeLimit 85 if config.strategy == 'GOA' or config.add_no_good_cuts: 86 config.logger.warning('ValueError: Cannot load a SolverResults object with bad status: error. ' 87 'MIP solver failed. This usually happens in the single-tree GOA algorithm. ' 88 "No-good cuts are added and GOA algorithm doesn't converge within the time limit. " 89 'No integer solution is found, so the cplex solver will report an error status. ') 90 return None, None 91 if main_mip_results.solver.termination_condition is tc.optimal: 92 if config.single_tree and not config.add_no_good_cuts and not regularization_problem: 93 if solve_data.objective_sense == minimize: 94 solve_data.LB = max( 95 main_mip_results.problem.lower_bound, solve_data.LB) 96 solve_data.bound_improved = solve_data.LB > solve_data.LB_progress[-1] 97 solve_data.LB_progress.append(solve_data.LB) 98 else: 99 solve_data.UB = min( 100 main_mip_results.problem.upper_bound, solve_data.UB) 101 solve_data.bound_improved = solve_data.UB < solve_data.UB_progress[-1] 102 solve_data.UB_progress.append(solve_data.UB) 103 104 elif main_mip_results.solver.termination_condition is tc.infeasibleOrUnbounded: 105 # Linear solvers will sometimes tell me that it's infeasible or 106 # unbounded during presolve, but fails to distinguish. We need to 107 # resolve with a solver option flag on. 108 main_mip_results, _ = distinguish_mip_infeasible_or_unbounded( 109 solve_data.mip, config) 110 return solve_data.mip, main_mip_results 111 112 if regularization_problem: 113 solve_data.mip.MindtPy_utils.objective_constr.deactivate() 114 solve_data.mip.MindtPy_utils.del_component('loa_proj_mip_obj') 115 solve_data.mip.MindtPy_utils.cuts.del_component('obj_reg_estimate') 116 if config.add_regularization == 'level_L1': 117 solve_data.mip.MindtPy_utils.del_component('L1_obj') 118 elif config.add_regularization == 'level_L_infinity': 119 solve_data.mip.MindtPy_utils.del_component( 120 'L_infinity_obj') 121 122 return solve_data.mip, main_mip_results 123 124 125def setup_mip_solver(solve_data, config, regularization_problem): 126 # Deactivate extraneous IMPORT/EXPORT suffixes 127 if config.nlp_solver == 'ipopt': 128 getattr(solve_data.mip, 'ipopt_zL_out', _DoNothing()).deactivate() 129 getattr(solve_data.mip, 'ipopt_zU_out', _DoNothing()).deactivate() 130 if regularization_problem: 131 mainopt = SolverFactory(config.mip_regularization_solver) 132 else: 133 if config.mip_solver == 'gurobi_persistent' and config.single_tree: 134 mainopt = GurobiPersistent4MindtPy() 135 mainopt.solve_data = solve_data 136 mainopt.config = config 137 else: 138 mainopt = SolverFactory(config.mip_solver) 139 140 # determine if persistent solver is called. 141 if isinstance(mainopt, PersistentSolver): 142 mainopt.set_instance(solve_data.mip, symbolic_solver_labels=True) 143 if config.single_tree and not regularization_problem: 144 # Configuration of cplex lazy callback 145 if config.mip_solver == 'cplex_persistent': 146 lazyoa = mainopt._solver_model.register_callback( 147 single_tree.LazyOACallback_cplex) 148 # pass necessary data and parameters to lazyoa 149 lazyoa.main_mip = solve_data.mip 150 lazyoa.solve_data = solve_data 151 lazyoa.config = config 152 lazyoa.opt = mainopt 153 mainopt._solver_model.set_warning_stream(None) 154 mainopt._solver_model.set_log_stream(None) 155 mainopt._solver_model.set_error_stream(None) 156 if config.mip_solver == 'gurobi_persistent': 157 mainopt.set_callback(single_tree.LazyOACallback_gurobi) 158 if config.use_tabu_list: 159 tabulist = mainopt._solver_model.register_callback( 160 tabu_list.IncumbentCallback_cplex) 161 tabulist.solve_data = solve_data 162 tabulist.opt = mainopt 163 tabulist.config = config 164 mainopt._solver_model.parameters.preprocessing.reduce.set(1) 165 # If the callback is used to reject incumbents, the user must set the 166 # parameter c.parameters.preprocessing.reduce either to the value 1 (one) 167 # to restrict presolve to primal reductions only or to 0 (zero) to disable all presolve reductions 168 mainopt._solver_model.set_warning_stream(None) 169 mainopt._solver_model.set_log_stream(None) 170 mainopt._solver_model.set_error_stream(None) 171 return mainopt 172 173 174# The following functions deal with handling the solution we get from the above MIP solver function 175 176 177def handle_main_optimal(main_mip, solve_data, config, update_bound=True): 178 """ 179 This function copies the result from 'solve_main' to the working model and updates the upper/lower bound. This 180 function is called after an optimal solution is found for the main problem. 181 182 Parameters 183 ---------- 184 main_mip: Pyomo model 185 the MIP main problem 186 solve_data: MindtPy Data Container 187 data container that holds solve-instance data 188 config: ConfigBlock 189 contains the specific configurations for the algorithm 190 """ 191 # proceed. Just need integer values 192 MindtPy = main_mip.MindtPy_utils 193 # check if the value of binary variable is valid 194 for var in MindtPy.discrete_variable_list: 195 if var.value is None: 196 config.logger.warning( 197 "Integer variable {} not initialized. It is set to it's lower bound".format(var.name)) 198 var.value = var.lb # nlp_var.bounds[0] 199 # warm start for the nlp subproblem 200 copy_var_list_values( 201 main_mip.MindtPy_utils.variable_list, 202 solve_data.working_model.MindtPy_utils.variable_list, 203 config) 204 205 if update_bound: 206 if solve_data.objective_sense == minimize: 207 solve_data.LB = max( 208 value(MindtPy.mip_obj.expr), solve_data.LB) 209 solve_data.bound_improved = solve_data.LB > solve_data.LB_progress[-1] 210 solve_data.LB_progress.append(solve_data.LB) 211 else: 212 solve_data.UB = min( 213 value(MindtPy.mip_obj.expr), solve_data.UB) 214 solve_data.bound_improved = solve_data.UB < solve_data.UB_progress[-1] 215 solve_data.UB_progress.append(solve_data.UB) 216 config.logger.info( 217 'MIP %s: OBJ: %s LB: %s UB: %s TIME: %ss' 218 % (solve_data.mip_iter, value(MindtPy.mip_obj.expr), 219 solve_data.LB, solve_data.UB, round(get_main_elapsed_time(solve_data.timing), 2))) 220 221 222def handle_main_other_conditions(main_mip, main_mip_results, solve_data, config): 223 """ 224 This function handles the result of the latest iteration of solving the MIP problem (given any of a few 225 edge conditions, such as if the solution is neither infeasible nor optimal). 226 227 Parameters 228 ---------- 229 main_mip: Pyomo model 230 the MIP main problem 231 main_mip_results: Pyomo results object 232 result from solving the MIP problem 233 solve_data: MindtPy Data Container 234 data container that holds solve-instance data 235 config: ConfigBlock 236 contains the specific configurations for the algorithm 237 """ 238 if main_mip_results.solver.termination_condition is tc.infeasible: 239 handle_main_infeasible(main_mip, solve_data, config) 240 elif main_mip_results.solver.termination_condition is tc.unbounded: 241 temp_results = handle_main_unbounded(main_mip, solve_data, config) 242 elif main_mip_results.solver.termination_condition is tc.infeasibleOrUnbounded: 243 temp_results = handle_main_unbounded(main_mip, solve_data, config) 244 if temp_results.solver.termination_condition is tc.infeasible: 245 handle_main_infeasible(main_mip, solve_data, config) 246 elif main_mip_results.solver.termination_condition is tc.maxTimeLimit: 247 handle_main_max_timelimit( 248 main_mip, main_mip_results, solve_data, config) 249 solve_data.results.solver.termination_condition = tc.maxTimeLimit 250 elif (main_mip_results.solver.termination_condition is tc.other and 251 main_mip_results.solution.status is SolutionStatus.feasible): 252 # load the solution and suppress the warning message by setting 253 # solver status to ok. 254 MindtPy = main_mip.MindtPy_utils 255 config.logger.info( 256 'MILP solver reported feasible solution, ' 257 'but not guaranteed to be optimal.') 258 copy_var_list_values( 259 main_mip.MindtPy_utils.variable_list, 260 solve_data.working_model.MindtPy_utils.variable_list, 261 config) 262 if solve_data.objective_sense == minimize: 263 solve_data.LB = max( 264 main_mip_results.problem.lower_bound, solve_data.LB) 265 solve_data.bound_improved = solve_data.LB > solve_data.LB_progress[-1] 266 solve_data.LB_progress.append(solve_data.LB) 267 else: 268 solve_data.UB = min( 269 main_mip_results.problem.upper_bound, solve_data.UB) 270 solve_data.bound_improved = solve_data.UB < solve_data.UB_progress[-1] 271 solve_data.UB_progress.append(solve_data.UB) 272 config.logger.info( 273 'MIP %s: OBJ: %s LB: %s UB: %s' 274 % (solve_data.mip_iter, value(MindtPy.mip_obj.expr), 275 solve_data.LB, solve_data.UB)) 276 else: 277 raise ValueError( 278 'MindtPy unable to handle MILP main termination condition ' 279 'of %s. Solver message: %s' % 280 (main_mip_results.solver.termination_condition, main_mip_results.solver.message)) 281 282 283def handle_main_infeasible(main_mip, solve_data, config): 284 """ 285 This function handles the result of the latest iteration of solving the MIP problem given an infeasible solution. 286 287 Parameters 288 ---------- 289 main_mip: Pyomo model 290 the MIP main problem 291 solve_data: MindtPy Data Container 292 data container that holds solve-instance data 293 config: ConfigBlock 294 contains the specific configurations for the algorithm 295 """ 296 config.logger.info( 297 'MILP main problem is infeasible. ' 298 'Problem may have no more feasible ' 299 'binary configurations.') 300 if solve_data.mip_iter == 1: 301 config.logger.warning( 302 'MindtPy initialization may have generated poor ' 303 'quality cuts.') 304 # TODO no-good cuts for single tree case 305 # set optimistic bound to infinity 306 if solve_data.objective_sense == minimize: 307 solve_data.LB_progress.append(solve_data.LB) 308 else: 309 solve_data.UB_progress.append(solve_data.UB) 310 config.logger.info( 311 'MindtPy exiting due to MILP main problem infeasibility.') 312 if solve_data.results.solver.termination_condition is None: 313 if solve_data.mip_iter == 0: 314 solve_data.results.solver.termination_condition = tc.infeasible 315 else: 316 solve_data.results.solver.termination_condition = tc.feasible 317 318 319def handle_main_max_timelimit(main_mip, main_mip_results, solve_data, config): 320 """ 321 This function handles the result of the latest iteration of solving the MIP problem given that solving the 322 MIP takes too long. 323 324 Parameters 325 ---------- 326 main_mip: Pyomo model 327 the MIP main problem 328 solve_data: MindtPy Data Container 329 data container that holds solve-instance data 330 config: ConfigBlock 331 contains the specific configurations for the algorithm 332 """ 333 # TODO check that status is actually ok and everything is feasible 334 MindtPy = main_mip.MindtPy_utils 335 config.logger.info( 336 'Unable to optimize MILP main problem ' 337 'within time limit. ' 338 'Using current solver feasible solution.') 339 copy_var_list_values( 340 main_mip.MindtPy_utils.variable_list, 341 solve_data.working_model.MindtPy_utils.variable_list, 342 config) 343 if solve_data.objective_sense == minimize: 344 solve_data.LB = max( 345 main_mip_results.problem.lower_bound, solve_data.LB) 346 solve_data.bound_improved = solve_data.LB > solve_data.LB_progress[-1] 347 solve_data.LB_progress.append(solve_data.LB) 348 else: 349 solve_data.UB = min( 350 main_mip_results.problem.upper_bound, solve_data.UB) 351 solve_data.bound_improved = solve_data.UB < solve_data.UB_progress[-1] 352 solve_data.UB_progress.append(solve_data.UB) 353 config.logger.info( 354 'MIP %s: OBJ: %s LB: %s UB: %s' 355 % (solve_data.mip_iter, value(MindtPy.mip_obj.expr), 356 solve_data.LB, solve_data.UB)) 357 358 359def handle_main_unbounded(main_mip, solve_data, config): 360 """ 361 This function handles the result of the latest iteration of solving the MIP problem given an unbounded solution 362 due to the relaxation. 363 364 Parameters 365 ---------- 366 main_mip: Pyomo model 367 the MIP main problem 368 solve_data: MindtPy Data Container 369 data container that holds solve-instance data 370 config: ConfigBlock 371 contains the specific configurations for the algorithm 372 """ 373 # Solution is unbounded. Add an arbitrary bound to the objective and resolve. 374 # This occurs when the objective is nonlinear. The nonlinear objective is moved 375 # to the constraints, and deactivated for the linear main problem. 376 MindtPy = main_mip.MindtPy_utils 377 config.logger.warning( 378 'main MILP was unbounded. ' 379 'Resolving with arbitrary bound values of (-{0:.10g}, {0:.10g}) on the objective. ' 380 'You can change this bound with the option obj_bound.'.format(config.obj_bound)) 381 MindtPy.objective_bound = Constraint( 382 expr=(-config.obj_bound, MindtPy.mip_obj.expr, config.obj_bound)) 383 mainopt = SolverFactory(config.mip_solver) 384 if isinstance(mainopt, PersistentSolver): 385 mainopt.set_instance(main_mip) 386 set_solver_options(mainopt, solve_data, config, solver_type='mip') 387 with SuppressInfeasibleWarning(): 388 main_mip_results = mainopt.solve( 389 main_mip, tee=config.mip_solver_tee, **config.mip_solver_args) 390 return main_mip_results 391 392 393def handle_regularization_main_tc(main_mip, main_mip_results, solve_data, config): 394 if main_mip_results is None: 395 config.logger.info( 396 'Failed to solve the regularization problem.' 397 'The solution of the OA main problem will be adopted.') 398 elif main_mip_results.solver.termination_condition in {tc.optimal, tc.feasible}: 399 handle_main_optimal( 400 main_mip, solve_data, config, update_bound=False) 401 elif main_mip_results.solver.termination_condition is tc.maxTimeLimit: 402 config.logger.info( 403 'Regularization problem failed to converge within the time limit.') 404 solve_data.results.solver.termination_condition = tc.maxTimeLimit 405 # break 406 elif main_mip_results.solver.termination_condition is tc.infeasible: 407 config.logger.info( 408 'Regularization problem infeasible.') 409 elif main_mip_results.solver.termination_condition is tc.unbounded: 410 config.logger.info( 411 'Regularization problem ubounded.' 412 'Sometimes solving MIQP in cplex, unbounded means infeasible.') 413 elif main_mip_results.solver.termination_condition is tc.unknown: 414 config.logger.info( 415 'Termination condition of the regularization problem is unknown.') 416 if main_mip_results.problem.lower_bound != float('-inf'): 417 config.logger.info('Solution limit has been reached.') 418 handle_main_optimal( 419 main_mip, solve_data, config, update_bound=False) 420 else: 421 config.logger.info('No solution obtained from the regularization subproblem.' 422 'Please set mip_solver_tee to True for more informations.' 423 'The solution of the OA main problem will be adopted.') 424 else: 425 raise ValueError( 426 'MindtPy unable to handle regularization problem termination condition ' 427 'of %s. Solver message: %s' % 428 (main_mip_results.solver.termination_condition, main_mip_results.solver.message)) 429 430 431def setup_main(solve_data, config, fp, regularization_problem): 432 """ 433 Set up main problem/main regularization problem for OA, ECP, Feasibility Pump and ROA methods. 434 435 Parameters 436 ---------- 437 solve_data: MindtPy Data Container 438 data container that holds solve-instance data 439 config: ConfigBlock 440 contains the specific configurations for the algorithm 441 fp: Bool 442 generate the feasibility pump regularization main problem 443 regularization_problem: Bool 444 generate the ROA regularization main problem 445 """ 446 MindtPy = solve_data.mip.MindtPy_utils 447 448 for c in MindtPy.constraint_list: 449 if c.body.polynomial_degree() not in {1, 0}: 450 c.deactivate() 451 452 MindtPy.cuts.activate() 453 454 sign_adjust = 1 if solve_data.objective_sense == minimize else - 1 455 MindtPy.del_component('mip_obj') 456 if regularization_problem and config.single_tree: 457 MindtPy.del_component('loa_proj_mip_obj') 458 MindtPy.cuts.del_component('obj_reg_estimate') 459 if config.add_regularization is not None and config.add_no_good_cuts: 460 if regularization_problem: 461 MindtPy.cuts.no_good_cuts.activate() 462 else: 463 MindtPy.cuts.no_good_cuts.deactivate() 464 465 if fp: 466 MindtPy.del_component('fp_mip_obj') 467 if config.fp_main_norm == 'L1': 468 MindtPy.fp_mip_obj = generate_norm1_objective_function( 469 solve_data.mip, 470 solve_data.working_model, 471 discrete_only=config.fp_discrete_only) 472 elif config.fp_main_norm == 'L2': 473 MindtPy.fp_mip_obj = generate_norm2sq_objective_function( 474 solve_data.mip, 475 solve_data.working_model, 476 discrete_only=config.fp_discrete_only) 477 elif config.fp_main_norm == 'L_infinity': 478 MindtPy.fp_mip_obj = generate_norm_inf_objective_function( 479 solve_data.mip, 480 solve_data.working_model, 481 discrete_only=config.fp_discrete_only) 482 elif regularization_problem: 483 if MindtPy.objective_list[0].expr.polynomial_degree() in {1, 0}: 484 MindtPy.objective_constr.activate() 485 if config.add_regularization == 'level_L1': 486 MindtPy.loa_proj_mip_obj = generate_norm1_objective_function(solve_data.mip, 487 solve_data.best_solution_found, 488 discrete_only=False) 489 elif config.add_regularization == 'level_L2': 490 MindtPy.loa_proj_mip_obj = generate_norm2sq_objective_function(solve_data.mip, 491 solve_data.best_solution_found, 492 discrete_only=False) 493 elif config.add_regularization == 'level_L_infinity': 494 MindtPy.loa_proj_mip_obj = generate_norm_inf_objective_function(solve_data.mip, 495 solve_data.best_solution_found, 496 discrete_only=False) 497 elif config.add_regularization in {'grad_lag', 'hess_lag', 'hess_only_lag', 'sqp_lag'}: 498 MindtPy.loa_proj_mip_obj = generate_lag_objective_function(solve_data.mip, 499 solve_data.best_solution_found, 500 config, 501 solve_data, 502 discrete_only=False) 503 if solve_data.objective_sense == minimize: 504 MindtPy.cuts.obj_reg_estimate = Constraint( 505 expr=MindtPy.objective_value <= (1 - config.level_coef) * solve_data.UB + config.level_coef * solve_data.LB) 506 else: 507 MindtPy.cuts.obj_reg_estimate = Constraint( 508 expr=MindtPy.objective_value >= (1 - config.level_coef) * solve_data.LB + config.level_coef * solve_data.UB) 509 510 else: 511 if config.add_slack: 512 MindtPy.del_component('aug_penalty_expr') 513 514 MindtPy.aug_penalty_expr = Expression( 515 expr=sign_adjust * config.OA_penalty_factor * sum( 516 v for v in MindtPy.cuts.slack_vars[...])) 517 main_objective = MindtPy.objective_list[-1] 518 MindtPy.mip_obj = Objective( 519 expr=main_objective.expr + 520 (MindtPy.aug_penalty_expr if config.add_slack else 0), 521 sense=solve_data.objective_sense) 522 523 if config.use_dual_bound: 524 # Delete previously added dual bound constraint 525 MindtPy.cuts.del_component('dual_bound') 526 if solve_data.objective_sense == minimize: 527 MindtPy.cuts.dual_bound = Constraint( 528 expr=main_objective.expr + 529 (MindtPy.aug_penalty_expr if config.add_slack else 0) >= solve_data.LB, 530 doc='Objective function expression should improve on the best found dual bound') 531 else: 532 MindtPy.cuts.dual_bound = Constraint( 533 expr=main_objective.expr + 534 (MindtPy.aug_penalty_expr if config.add_slack else 0) <= solve_data.UB, 535 doc='Objective function expression should improve on the best found dual bound') 536