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# NOTE: this solver is disabled (see the first try block below). This 12# code is out of date, and this code is not regularly tested by Pyomo developers. 13# The python-glpk package is only supported on Debian Linux platforms, so 14# it is not clear if this is a valuable solver interface, particularly since 15# commercial vendors now have good support for Python interfaces (e.g. CPLEX and 16# Gurobi). 17 18import sys 19_glpk_version = None 20glpk_python_api_exists = None 21def configure_glpk_direct(): 22 global _glpk_version 23 global glpk_python_api_exists 24 try: 25 # import all the glp_* functions 26 if False: # DISABLED 27 import glpk 28 glpk_python_api_exists = True 29 else: 30 glpk_python_api_exists = False 31 except ImportError: 32 glpk_python_api_exists = False 33 except Exception as e: 34 # other forms of exceptions can be thrown by the glpk python 35 # import. For example, an error in code invoked by the module's 36 # __init__. We should continue gracefully and not cause a fatal 37 # error in Pyomo. 38 print("Import of glpk failed - glpk message="+str(e)+"\n") 39 glpk_python_api_exists = False 40 41from pyomo.common.collections import Bunch 42 43from pyomo.opt.base import OptSolver 44from pyomo.opt.base.solvers import _extract_version, SolverFactory 45from pyomo.opt.results import SolverResults, SolverStatus, SolutionStatus, ProblemSense 46from pyomo.core.base.numvalue import value 47 48import logging 49logger = logging.getLogger('pyomo.solvers') 50 51 52 53@SolverFactory.register('_glpk_direct', doc='Direct Python interface to the GLPK LP/MIP solver.') 54class GLPKDirect ( OptSolver ): 55 """The GLPK LP/MIP solver (direct API plugin) 56 57 The glpk_direct plugin offers an API interface to the GLPK. It requires the 58 Python-GLPK API interface provided via the SWIG interface available through 59 most Linux distributions repositories. For more information, see Jo~ao Pedro 60 Pedroso's (the author) page at http://www.dcc.fc.up.pt/~jpp/ 61 62 Because of the direct connection with the GLPK API, no temporary files need be 63 written or read. That ostensibly makes this a faster plugin than the 64 file-based glpk plugin. However, you will likely not notice any speed up 65 unless you are using the GLPK solver with PySP problems (due to the rapid 66 re-solves). 67 68 One downside to the lack of temporary files, is that there is no LP file to 69 inspect for clues while debugging a model. For that, use the write_lp solver 70 option: 71 72 $ pyomo model.{py,dat} \ 73 --solver=glpk_direct \ 74 --solver-options write_lp=/path/to/some/file.lp 75 76 One can also specify the particular GLPK algorithm to use with the 'algorithm' 77 solver option. There are 4 algorithms: 78 79 simplex - the default algorithm for non-MIP problems (primal simplex) 80 intopt - the default algorithm for MIP problems 81 exact - tentative implementation of two-phase primal simplex with exact 82 arithmetic internally. 83 interior - only aplicable to LP problems 84 85 $ pyomo model.{py,dat} \ 86 --solver glpk_direct \ 87 --solver-options algorithm=exact 88 89 For more information on available algorithms, see the GLPK documentation. 90 """ 91 92 def __init__(self, **kwds): 93 # 94 # Call base class constructor 95 # 96 kwds['type'] = 'glpk_direct' 97 OptSolver.__init__(self, **kwds) 98 99 # NOTE: eventually both of the following attributes should be migrated 100 # to a common base class. Is the current solve warm-started? A 101 # transient data member to communicate state information across the 102 # _presolve, _apply_solver, and _postsolve methods. 103 self.warm_start_solve = False 104 self._timelimit = None 105 106 # Note: Undefined capabilities default to 'None' 107 self._capabilities = Bunch() 108 self._capabilities.linear = True 109 self._capabilities.integer = True 110 111 def _get_version(self): 112 """ 113 Returns a tuple describing the solver executable version. 114 """ 115 if _glpk_version is None: 116 return _extract_version('') 117 return _glpk_version 118 119 def _populate_glpk_instance ( self, model ): 120 121 from pyomo.core.base import Var, Objective, Constraint, SOSConstraint 122 123 try: 124 lp = glp_create_prob() 125 except Exception: 126 e = sys.exc_info()[1] 127 msg = 'Unable to create GLPK problem instance. Have you installed' \ 128 '\n the Python bindings for GLPK?\n\n\tError message: %s' 129 raise Exception(msg % e) 130 131 objective = sorted( model.component_map(Objective, active=True).values() )[0] 132 # so we can correctly map the solution to the correct objective label in _postsolve 133 lp.objective_name = sorted( model.component_map(Objective, active=True).keys() )[0] 134 sense = GLP_MAX 135 if objective.is_minimizing(): sense = GLP_MIN 136 137 constraint_list = model.component_map(Constraint, active=True) 138 variable_list = model.component_map(Var) 139 num_constraints = model.statistics.number_of_constraints 140 num_variables = model.statistics.number_of_variables 141 142 sosn = self._capabilities.sosn 143 sos1 = self._capabilities.sos1 144 sos2 = self._capabilities.sos2 145 146 for soscondata in model.component_data_objects(SOSConstraint, active=True): 147 raise Exception("Solver: glpk_direct does not support SOSConstraint declarations") 148 149 glp_set_prob_name(lp, model.name) 150 151 glp_set_obj_dir( lp, sense ) 152 glp_add_rows( lp, num_constraints ) 153 glp_add_cols( lp, num_variables ) 154 155 # 1 extra because GLPK's arrays in this context are 1-based, not 0-based 156 coef_count = num_constraints * num_variables + 1 157 Ai = intArray( coef_count ) 158 Aj = intArray( coef_count ) 159 Ar = doubleArray( coef_count ) 160 161 row = col = coef_count = 0 162 colvar_map = dict() 163 rowvar_map = dict() 164 165 # In matrix parlance, variables are columns 166 for name in variable_list: 167 var_set = variable_list[ name ] 168 for ii in var_set: 169 var = var_set[ ii ] 170 if var.fixed is True: 171 continue 172 173 lb = ub = 0.0 174 if (not var.has_lb()) and \ 175 (not var.has_ub()): 176 var_type = GLP_FR 177 elif not var.has_lb(): 178 var_type = GLP_UB 179 ub = value(var.ub) 180 elif not var.has_ub(): 181 var_type = GLP_LO 182 lb = value(var.lb) 183 else: 184 var_type = GLP_DB 185 lb = value(var.lb) 186 ub = value(var.ub) 187 188 col += 1 189 colvar_map[ var.label ] = col 190 191 # the name is perhaps not necessary, but for completeness ... 192 glp_set_col_name( lp, col, var.label ) 193 glp_set_col_bnds( lp, col, var_type, lb, ub ) 194 195 # Be sure to impart the integer and binary nature of any variables 196 if var.is_integer(): 197 glp_set_col_kind( lp, col, GLP_IV ) 198 elif var.is_binary(): 199 glp_set_col_kind( lp, col, GLP_BV ) 200 elif var.is_continuous(): 201 glp_set_col_kind( lp, col, GLP_CV ) # continuous 202 else: 203 raise TypeError("Invalid domain type for variable with name '%s'. " 204 "Variable is not continuous, integer, or binary.") 205 206 model_canonical_repn = getattr(model, "_canonical_repn", None) 207 if model_canonical_repn is None: 208 raise ValueError("No _canonical_repn ComponentMap was found on " 209 "block with name %s. Did you forget to preprocess?" 210 % (model.name)) 211 212 for name in constraint_list: 213 constraint_set = constraint_list[ name ] 214 215 for ii in constraint_set: 216 constraint = constraint_set[ ii ] 217 if not constraint.active: continue 218 elif (not constraint.has_lb()) and \ 219 (not constraint.has_ub()): 220 continue 221 222 expression = model_canonical_repn.get(constraint) 223 if constraint is None: 224 raise ValueError("No entry found in _canonical_repn ComponentMap on " 225 "block %s for active constraint with name %s. " 226 "Did you forget to preprocess?" 227 % (model.name, constraint.name)) 228 229 offset = 0.0 230 if 0 in expression: 231 offset = expression[0][None] 232 233 lbound = ubound = -offset 234 235 if constraint.equality: 236 var_type = GLP_FX # Fixed 237 lbound = ubound = constraint.lower() - offset 238 elif not constraint.has_lb(): 239 var_type = GLP_UP # Upper bounded only 240 ubound += constraint.upper() 241 elif not constraint.has_ub(): 242 var_type = GLP_LO # Lower bounded only 243 lbound += constraint.lower() 244 else: 245 var_type = GLP_DB # Double bounded 246 lbound += constraint.lower() 247 ubound += constraint.upper() 248 249 row += 1 250 rowvar_map[ constraint.label ] = row 251 252 # just as with variables, set the name just for completeness ... 253 glp_set_row_name( lp, row, constraint.label ) 254 glp_set_row_bnds( lp, row, var_type, lbound, ubound ) 255 256 if 1 in expression: # first-order terms 257 keys = sorted( expression[1].keys() ) 258 for var_key in keys: 259 index = var_key.keys()[0] 260 var = expression[-1][ index ] 261 coef = expression[ 1][ var_key ] 262 col = colvar_map[ var.label ] 263 264 coef_count += 1 265 Ai[ coef_count ] = row 266 Aj[ coef_count ] = col 267 Ar[ coef_count ] = coef 268 269 # with the rows and columns named and bounded, load the coefficients 270 glp_load_matrix( lp, coef_count, Ai, Aj, Ar ) 271 272 for key in objective: 273 274 expression = model_canonical_repn.get(objective[key]) 275 if expression is None: 276 raise ValueError("No entry found in _canonical_repn ComponentMap on " 277 "block %s for active objective with name %s. " 278 "Did you forget to preprocess?" 279 % (model.name, objective[key].name)) 280 281 if expression.is_constant(): 282 msg = "Ignoring objective '%s[%s]' which is constant" 283 logger.warning( msg % (str(objective), str(key)) ) 284 continue 285 286 if 1 in expression: # first-order terms 287 keys = sorted( expression[1].keys() ) 288 for var_key in keys: 289 index = var_key.keys()[0] 290 label = expression[-1][ index ].label 291 coef = expression[ 1][ var_key ] 292 col = colvar_map[ label ] 293 glp_set_obj_coef( lp, col, coef ) 294 295 elif -1 in expression: 296 pass 297 else: 298 msg = "Nonlinear objective to GLPK. GLPK can only handle " \ 299 "linear problems." 300 raise RuntimeError( msg ) 301 302 303 self._glpk_instance = lp 304 self._glpk_rowvar_map = rowvar_map 305 self._glpk_colvar_map = colvar_map 306 307 308 def warm_start_capable(self): 309 # Note to dev who gets back here: GLPK has the notion of presolving, but 310 # it's "built-in" to each optimization function. To disable it, make use 311 # of the second argument of glp_smcp type. See the GLPK documentation 312 # PDF for further information. 313 msg = "GLPK has the ability to use warmstart solutions. However, it " \ 314 "has not yet been implemented into the Pyomo glpk_direct plugin." 315 logger.info( msg ) 316 return False 317 318 319 def warm_start(self, instance): 320 pass 321 322 323 def _presolve(self, *args, **kwargs): 324 from pyomo.core.base.PyomoModel import Model 325 326 self.warm_start_solve = kwargs.pop( 'warmstart', False ) 327 328 model = args[0] 329 if len(args) != 1: 330 msg = "The glpk_direct plugin method '_presolve' must be supplied " \ 331 "a single problem instance - %s were supplied" 332 raise ValueError(msg % len(args)) 333 elif not isinstance(model, Model): 334 msg = "The problem instance supplied to the glpk_direct plugin " \ 335 "'_presolve' method must be of type 'Model'" 336 raise ValueError(msg) 337 338 self._populate_glpk_instance( model ) 339 lp = self._glpk_instance 340 self.is_integer = ( glp_get_num_int( lp ) > 0 and True or False ) 341 342 if 'write_lp' in self.options: 343 fname = self.options.write_lp 344 glp_write_lp( lp, None, fname ) 345 346 self.algo = 'Simplex (primal)' 347 algorithm = glp_simplex 348 if self.is_integer > 0: 349 self.algo = 'Mixed Integer' 350 algorithm = glp_intopt 351 352 if 'algorithm' in self.options: 353 if 'simplex' == self.options.algorithm: 354 self.algo = 'Simplex (primal)' 355 algorithm = glp_simplex 356 elif 'exact' == self.options.algorithm: 357 self.algo = 'Simplex (two-phase primal)' 358 algorithm = glp_exact 359 elif 'interior' == self.options.algorithm: 360 self.algo = 'Interior Point' 361 algorithm = glp_interior 362 elif 'intopt' == self.options.algorithm: 363 self.alpo = 'Mixed Integer' 364 algorithm = glp_intopt 365 else: 366 msg = "Unknown solver specified\n Unknown: %s\n Using: %s\n" 367 logger.warning( msg % (self.options.algorithm, self.algo) ) 368 self._algorithm = algorithm 369 370 if 'Simplex (primal)' == self.algo: 371 parm = glp_smcp() 372 glp_init_smcp( parm ) 373 elif 'Simplex (two-phase primal)' == self.algo: 374 parm = glp_smcp() 375 glp_init_smcp( parm ) 376 elif 'Interior Point' == self.algo: 377 parm = glp_iptcp() 378 glp_init_iptcp( parm ) 379 elif 'Mixed Integer' == self.algo: 380 parm = glp_iocp() 381 glp_init_iocp( parm ) 382 if self.options.mipgap: 383 parm.mip_gap = self.options.mipgap 384 385 if self._timelimit and self._timelimit > 0.0: 386 parm.tm_lim = self._timelimit 387 388 parm.msg_lev = GLP_MSG_OFF 389 parm.presolve = GLP_ON 390 391 self._solver_params = parm 392 393 # Scaffolding in place: I /believe/ GLPK can do warmstarts; just supply 394 # the basis solution and turn of the presolver 395 if self.warm_start_solve is True: 396 397 if len(args) != 1: 398 msg = "The glpk_direct _presolve method can only handle a single "\ 399 "problem instance - %s were supplied" 400 raise ValueError(msg % len(args)) 401 402 parm.presolve = GLP_OFF 403 self.warm_start( model ) 404 405 406 def _apply_solver(self): 407 lp = self._glpk_instance 408 parm = self._solver_params 409 algorithm = self._algorithm 410 411 # Actually solve the problem. 412 try: 413 beg = glp_time() 414 self.solve_return_code = algorithm( self._glpk_instance, parm ) 415 end = glp_time() 416 self._glpk_solve_time = glp_difftime( end, beg ) 417 except Exception: 418 e = sys.exc_info()[1] 419 msg = str(e) 420 if 'algorithm' in self.options: 421 msg = "Unexpected error using '%s' algorithm. Is it the correct "\ 422 "correct algorithm for the problem type?" 423 msg %= self.options.algorithm 424 logger.error( msg ) 425 raise 426 427 # FIXME: can we get a return code indicating if GLPK had a 428 # significant failure? 429 return Bunch(rc=None, log=None) 430 431 432 def _glpk_return_code_to_message ( self ): 433 code = self.solve_return_code 434 if 0 == code: 435 return "Algorithm completed successfully. (This does not " \ 436 "necessarily mean an optimal solution was found.)" 437 elif GLP_EBADB == code: 438 return "Unable to start the search, because the initial basis " \ 439 "specified in the problem object is invalid -- the number of " \ 440 "basic (auxiliary and structural) variables is not the same as " \ 441 "the number of rows in the problem object." 442 elif GLP_ESING == code: 443 return "Unable to start the search, because the basis matrix " \ 444 "corresponding to the initial basis is singular within the " \ 445 "working precision." 446 elif GLP_ECOND == code: 447 return "Unable to start the search, because the basis matrix " \ 448 "corresponding to the initial basis is ill-conditioned, i.e. its "\ 449 "condition number is too large." 450 elif GLP_EBOUND == code: 451 return "Unable to start the search, because some double-bounded " \ 452 "(auxiliary or structural) variables have incorrect bounds." 453 elif GLP_EFAIL == code: 454 return "The search was prematurely terminated due to the solver " \ 455 "failure." 456 elif GLP_EOBJLL == code: 457 return "The search was prematurely terminated, because the " \ 458 "objective function being maximized has reached its lower limit " \ 459 "and continues decreasing (the dual simplex only)." 460 elif GLP_EOBJUL == code: 461 return "The search was prematurely terminated, because the " \ 462 "objective function being minimized has reached its upper limit " \ 463 "and continues increasing (the dual simplex only)." 464 elif GLP_EITLIM == code: 465 return "The search was prematurely terminated, because the simplex " \ 466 "iteration limit has been exceeded." 467 elif GLP_ETMLIM == code: 468 return "The search was prematurely terminated, because the time " \ 469 "limit has been exceeded." 470 elif GLP_ENOPFS == code: 471 return "The LP problem instance has no primal feasible solution " \ 472 "(only if the LP presolver is used)." 473 elif GLP_ENODFS == code: 474 return "The LP problem instance has no dual feasible solution " \ 475 "(only if the LP presolver is used)." 476 else: 477 return "Unexpected error condition. Please consider remitting " \ 478 "this problem to the Pyomo developers and/or the GLPK project " \ 479 "so they can improve their softwares." 480 481 482 def _glpk_get_solution_status ( self ): 483 getstatus = glp_get_status 484 if self.is_integer: getstatus = glp_mip_status 485 486 status = getstatus( self._glpk_instance ) 487 if GLP_OPT == status: return SolutionStatus.optimal 488 elif GLP_FEAS == status: return SolutionStatus.feasible 489 elif GLP_INFEAS == status: return SolutionStatus.infeasible 490 elif GLP_NOFEAS == status: return SolutionStatus.other 491 elif GLP_UNBND == status: return SolutionStatus.other 492 elif GLP_UNDEF == status: return SolutionStatus.other 493 raise RuntimeError("Unknown solution status returned by GLPK solver") 494 495 496 def _glpk_get_solver_status ( self ): 497 rc = self.solve_return_code 498 if 0 == rc: return SolverStatus.ok 499 elif GLP_EBADB == rc: return SolverStatus.error 500 elif GLP_ESING == rc: return SolverStatus.error 501 elif GLP_ECOND == rc: return SolverStatus.error 502 elif GLP_EBOUND == rc: return SolverStatus.error 503 elif GLP_EFAIL == rc: return SolverStatus.aborted 504 elif GLP_EOBJLL == rc: return SolverStatus.aborted 505 elif GLP_EOBJUL == rc: return SolverStatus.aborted 506 elif GLP_EITLIM == rc: return SolverStatus.aborted 507 elif GLP_ETMLIM == rc: return SolverStatus.aborted 508 elif GLP_ENOPFS == rc: return SolverStatus.warning 509 elif GLP_ENODFS == rc: return SolverStatus.warning 510 else: return SolverStatus.unkown 511 512 513 def _postsolve(self): 514 lp = self._glpk_instance 515 num_variables = glp_get_num_cols( lp ) 516 bin_variables = glp_get_num_bin( lp ) 517 int_variables = glp_get_num_int( lp ) 518 519 # check suffixes 520 for suffix in self._suffixes: 521 if True: 522 raise RuntimeError("***The glpk_direct solver plugin cannot extract solution suffix="+suffix) 523 524 525 tpeak = glp_long() 526 glp_mem_usage( None, None, None, tpeak ) 527 # black magic trickery, thanks to Python's lack of pointers and SWIG's 528 # automatic API conversion 529 peak_mem = tpeak.lo 530 531 results = SolverResults() 532 soln = Solution() 533 prob = results.problem 534 solv = results.solver 535 536 solv.name = "GLPK " + glp_version() 537 solv.status = self._glpk_get_solver_status() 538 solv.return_code = self.solve_return_code 539 solv.message = self._glpk_return_code_to_message() 540 solv.algorithm = self.algo 541 solv.memory_used = "%d bytes, (%d KiB)" % (peak_mem, peak_mem/1024) 542 # solv.user_time = None 543 # solv.system_time = None 544 solv.wallclock_time = self._glpk_solve_time 545 # solv.termination_condition = None 546 # solv.termination_message = None 547 548 prob.name = glp_get_prob_name(lp) 549 prob.number_of_constraints = glp_get_num_rows(lp) 550 prob.number_of_nonzeros = glp_get_num_nz(lp) 551 prob.number_of_variables = num_variables 552 prob.number_of_binary_variables = bin_variables 553 prob.number_of_integer_variables = int_variables 554 prob.number_of_continuous_variables = num_variables - int_variables 555 prob.number_of_objectives = 1 556 557 prob.sense = ProblemSense.minimize 558 if GLP_MAX == glp_get_obj_dir( lp ): 559 prob.sense = ProblemSense.maximize 560 561 soln.status = self._glpk_get_solution_status() 562 563 if soln.status in ( SolutionStatus.optimal, SolutionStatus.feasible ): 564 get_col_prim = glp_get_col_prim 565 get_row_prim = glp_get_row_prim 566 get_obj_val = glp_get_obj_val 567 if self.is_integer: 568 get_col_prim = glp_mip_col_val 569 get_row_prim = glp_mip_row_val 570 get_obj_val = glp_mip_obj_val 571 572 obj_val = get_obj_val( lp ) 573 if prob.sense == ProblemSense.minimize: 574 prob.lower_bound = obj_val 575 else: 576 prob.upper_bound = obj_val 577 578 objective_name = lp.objective_name 579 soln.objective[objective_name] = {'Value': obj_val} 580 581 colvar_map = self._glpk_colvar_map 582 rowvar_map = self._glpk_rowvar_map 583 584 for var_label in colvar_map: 585 col = colvar_map[ var_label ] 586 soln.variable[ var_label ] = {"Value" : get_col_prim( lp, col )} 587 588 for row_label in rowvar_map: 589 row = rowvar_map[ row_label ] 590 soln.constraint[ row_label ] = {"Value" : get_row_prim( lp, row )} 591 592 results.solution.insert(soln) 593 594 self.results = results 595 596 # All done with the GLPK object, so free up some memory. 597 glp_free( lp ) 598 del self._glpk_instance, lp 599 600 # let the base class deal with returning results. 601 return OptSolver._postsolve(self) 602 603 604# TODO: add MockGLPKDirect class 605 606if not glpk_python_api_exists: 607 SolverFactory().unregister('_glpk_direct') 608 # SolverFactory().deactivate('_mock_glpk_direct') 609 610# vim: set fileencoding=utf-8 611