# _________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright (c) 2014 Sandia Corporation. # Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, # the U.S. Government retains certain rights in this software. # This software is distributed under the BSD License. # _________________________________________________________________________ # TODO: have ddsip convert create symbols files for second stage # TODO: parse second-stage solution and load into scenario tree workers # TODO: objective, cost, stage_costs # TODO: make output_scenario_tree_solution work import io import os import sys import time import shutil import logging import traceback import pyutilib.subprocess import pyutilib.services from pyomo.core import maximize from pyomo.opt import (TerminationCondition, SolverStatus, SolutionStatus) from pyomo.pysp.util.configured_object import PySPConfiguredObject from pyomo.pysp.util.config import (PySPConfigValue, PySPConfigBlock, safe_register_common_option, safe_register_unique_option, safe_declare_common_option, safe_declare_unique_option, _domain_tuple_of_str_or_dict, _domain_must_be_str) from pyomo.pysp.util.misc import (parse_command_line, launch_command) from pyomo.pysp.scenariotree.manager import \ (InvocationType, ScenarioTreeManagerFactory) import pyomo.pysp.convert.ddsip from pyomo.pysp.embeddedsp import EmbeddedSP from pyomo.pysp.solvers.spsolver import (SPSolverResults, SPSolverFactory) from pyomo.pysp.solvers.spsolvershellcommand import \ SPSolverShellCommand logger = logging.getLogger('pyomo.pysp') thisfile = os.path.abspath(__file__) _ddsip_group_label = "DDSIP Options" _firststage_var_suffix = "__DDSIP_FIRSTSTAGE" # maps the ddsip status code to tuples of the form # (SolutionStatus, SolverStatus, TerminationCondition, message) # feel free to modify if you have opinions one these _ddsip_status_map = {} # the process was terminated by the user _ddsip_status_map[-1] = (SolutionStatus.unknown, SolverStatus.aborted, TerminationCondition.userInterrupt, "Termination signal received.") # the node limit has been reached _ddsip_status_map[1] = (SolutionStatus.stoppedByLimit, SolverStatus.ok, TerminationCondition.maxEvaluations, "Node limit reached (total number of nodes).") # the gap (absolute or relative) has been reached _ddsip_status_map[2] = (SolutionStatus.optimal, SolverStatus.ok, TerminationCondition.optimal, "Gap reached.") # the time limit has been reached _ddsip_status_map[3] = (SolutionStatus.stoppedByLimit, SolverStatus.ok, TerminationCondition.maxTimeLimit, "Time limit reached.") # the maximal dispersion, i.e. the maximal difference of the # first-stage components within all remaining front nodes, # was less then the parameter NULLDISP (null dispersion) _ddsip_status_map[4] = (SolutionStatus.optimal, SolverStatus.ok, TerminationCondition.minStepLength, "Maximal dispersion equals zero.") # the whole branching tree was backtracked. _ddsip_status_map[5] = (SolutionStatus.optimal, SolverStatus.ok, TerminationCondition.minStepLength, ("The whole branching tree was " "backtracked. Probably due to MIP " "gaps (see below) the specified " "gap tolerance could not be reached.")) # no valid lower bound _ddsip_status_map[6] = (SolutionStatus.unknown, SolverStatus.warning, TerminationCondition.invalidProblem, "No valid lower bound.") # problem infeasible _ddsip_status_map[7] = (SolutionStatus.infeasible, SolverStatus.warning, TerminationCondition.infeasible, "Problem infeasible.") # problem unbounded _ddsip_status_map[8] = (SolutionStatus.unbounded, SolverStatus.warning, TerminationCondition.unbounded, "Problem unbounded.") def _load_solution(manager, scenario, solution_filename, info_filename, firststage_symbols_filename, secondstage_symbols_filename, scenario_id): # parse the symbol map firststage_symbol_map = {} with open(firststage_symbols_filename) as f: for line in f: lp_symbol, scenario_tree_id = line.strip().split() firststage_symbol_map[lp_symbol] = scenario_tree_id secondstage_symbol_map = {} with open(secondstage_symbols_filename) as f: for line in f: lp_symbol, scenario_tree_id = line.strip().split() secondstage_symbol_map[lp_symbol] = scenario_tree_id x_firststage = scenario._x[scenario.node_list[0].name] assert scenario.node_list[0] is manager.scenario_tree.findRootNode() x_secondstage = scenario._x[scenario.node_list[1].name] with open(solution_filename, 'r') as f: line = f.readline() while line.strip() != "1. Best Solution": line = f.readline() line = f.readline() assert line.startswith("Variable name Value") line = f.readline() assert line.startswith("-----------------------------------") line = f.readline().strip() while line != "": varlabel, varsol = line.split() x_firststage[firststage_symbol_map[varlabel]] = \ float(varsol) line = f.readline().strip() while line != "4. Second-stage solutions": line = f.readline().strip() scenario_label = ("Scenario %d:" % scenario_id) line = f.readline().strip() while line != scenario_label: line = f.readline().strip() line = f.readline().strip() while (line != "") and (not line.startswith("Scenario ")): varlabel, varsol = line.split() if varlabel != "ONE_VAR_CONSTANT": x_secondstage[secondstage_symbol_map[varlabel]] = \ float(varsol) line = f.readline().strip() return x class DDSIPSolver(SPSolverShellCommand, PySPConfiguredObject): @classmethod def _declare_options(cls, options=None): if options is None: options = PySPConfigBlock() safe_declare_unique_option( options, "firststage_suffix", PySPConfigValue( "__DDSIP_FIRSTSTAGE", domain=_domain_must_be_str, description=( "The suffix used to identity first-stage " "variables to DDSIP. Default is " "'__DDSIP_FIRSTSTAGE'" ), doc=None, visibility=0), ap_group=_ddsip_group_label) safe_declare_unique_option( options, "config_file", PySPConfigValue( None, domain=_domain_must_be_str, description=( "The name of a partial DDSIP configuration file " "that contains option specifications unrelated to " "the problem structure. If specified, the contents " "of this file will be appended to the " "configuration created by this solver interface. " "Default is None." ), doc=None, visibility=0), ap_group=_ddsip_group_label) safe_declare_common_option(options, "verbose", ap_group=_ddsip_group_label) return options def __init__(self): super(DDSIPSolver, self).__init__(self.register_options()) self.set_options_to_default() self._executable = "ddsip" def set_options_to_default(self): self._options = self.register_options() self._options._implicit_declaration = True @property def options(self): return self._options @property def name(self): return "ddsip" def _solve_impl(self, sp, output_solver_log=False, **kwds): """ Solve a stochastic program with the DDSIP solver. See the 'solve' method on the base class for additional keyword documentation. Args: sp: The stochastic program to solve. output_solver_log (bool): Stream the solver output during the solve. **kwds: Passed to the DDSIP file writer as I/O options (e.g., symbolic_solver_labels=True). Returns: A results object with information about the solution. """ if len(sp.scenario_tree.stages) > 2: raise ValueError("DDSIP solver does not handle more " "than 2 time-stages") # # Setup the DDSIP working directory # working_directory = self._create_tempdir("workdir") input_directory = os.path.join(working_directory, "ddsip_files") output_directory = os.path.join(input_directory, "sipout") logfile = self._files['logfile'] = \ os.path.join(working_directory, "ddsip.log") os.makedirs(input_directory) assert os.path.exists(input_directory) assert not os.path.exists(output_directory) info_filename = os.path.join(output_directory, "sip.out") solution_filename = os.path.join(output_directory, "solution.out") # # Create the DDSIP input files # if self.get_option("verbose"): print("Writing solver files in directory: %s" % (working_directory)) input_files = pyomo.pysp.convert.ddsip.\ convert_external( input_directory, self.options.firststage_suffix, sp, io_options=kwds) for key in input_files: self._add_tempfile(key, input_files[key]) self._update_config(input_files["config"]) # # Launch DDSIP # _cmd_string = self.executable+" < "+input_files["script"] if self.get_option("verbose"): print("Launching DDSIP solver with command: %s" % (_cmd_string)) ddsipstdin = None with open(input_files["script"]) as f: ddsipstdin = f.read() assert ddsipstdin is not None start = time.time() rc, log = pyutilib.subprocess.run( self.executable, cwd=input_directory, stdin=ddsipstdin, outfile=logfile, tee=output_solver_log) stop = time.time() if rc: if not self.available(): raise ValueError( "Solver executable does not exist: '%s'. " "(note that the default executable generated " "by the DDSIP build system does not have this name)" % (self.executable)) else: raise RuntimeError( "A nonzero return code (%s) was encountered after " "launching the command: %s. Check the solver log file " "for more information: %s" % (rc, _cmd_string, logfile)) # # Parse the DDSIP solution # if self.get_option("verbose"): print("Reading DDSIP solution from file: %s" % (solution_filename)) assert os.path.exists(output_directory) async_xhat = None async_responses = [] # TODO: The solution symbols are not being translated # back to scenario tree ids. Fix this and also have # the first and second stage solutions be loaded into # the Pyomo models #for scenario_id, scenario in enumerate(sp.scenario_tree.scenarios, 1): # async_responses.append(sp.invoke_function( # "_load_solution", # thisfile, # invocation_type=InvocationType.OnScenario(scenario.name), # function_args=(solution_filename, scenario_id), # async_call=True)) xhat, results = self._read_solution(sp, input_files["symbols"], info_filename, solution_filename) results.xhat = None if xhat is not None: results.xhat = {sp.scenario_tree.findRootNode().name: xhat} for res in async_responses: res.complete() return results def _update_config(self, config_filename): """ Writes a DDSIP config file """ # find the byte position where # we start appending to the config file # (just before END) append_pos = 0 with io.open(config_filename, mode='rb', buffering=0) as f: f.seek(0) append_pos = f.tell() for line in f: if line.strip().decode() == "END": break append_pos = f.tell() assert append_pos > 0 config_lines = {} config_lines[None] = [] config_lines['CPLEX'] = [] config_lines['CPLEXEEV'] = [] config_lines['CPLEXLB'] = [] config_lines['CPLEX2LB'] = [] config_lines['CPLEXUB'] = [] config_lines['CPLEX2UB'] = [] config_lines['CPLEXDUAL'] = [] config_lines['CPLEX2DUAL'] = [] has_cplex_opts = False # parse the user specified config file if self.options.config_file is not None: with open(self.options.config_file) as fconfig: section = config_lines[None] for line in fconfig: stripped = line.strip() if (stripped == "BEGIN") and \ (stripped == "END"): continue if "CPLEXBEGIN" in stripped: section = config_lines['CPLEX'] has_cplex_opts = True continue elif "CPLEXEEV" in stripped: section = config_lines['CPLEXEEV'] has_cplex_opts = True continue elif "CPLEXLB" in stripped: section = config_lines['CPLEXLB'] has_cplex_opts = True continue elif "CPLEX2LB" in stripped: section = config_lines['CPLEX2LB'] has_cplex_opts = True continue elif "CPLEXUB" in stripped: section = config_lines['CPLEXUB'] has_cplex_opts = True continue elif "CPLEX2UB" in stripped: section = config_lines['CPLEX2UB'] has_cplex_opts = True continue elif "CPLEXDUAL" in stripped: section = config_lines['CPLEXDUAL'] has_cplex_opts = True continue elif "CPLEX2DUAL" in stripped: section = config_lines['CPLEX2DUAL'] has_cplex_opts = True continue elif "CPLEXEND" in stripped: section = config_lines[None] continue section.append(line) for key in self.options: if (key != "firststage_suffix") and \ (key != "config_file") and \ (key != "verbose"): if key.startswith("CPLEX_"): section = config_lines["CPLEX"] prefix = "CPLEX_" has_cplex_opts = True elif key.startswith("CPLEXEEV_"): section = config_lines["CPLEXEEV"] prefix = "CPLEXEEV_" has_cplex_opts = True elif key.startswith("CPLEXLB_"): section = config_lines["CPLEXLB"] prefix = "CPLEXLB_" has_cplex_opts = True elif key.startswith("CPLEX2LB_"): section = config_lines["CPLEX2LB"] prefix = "CPLEX2LB_" has_cplex_opts = True elif key.startswith("CPLEXUB_"): section = config_lines["CPLEXUB"] prefix = "CPLEXUB_" has_cplex_opts = True elif key.startswith("CPLEX2UB_"): section = config_lines["CPLEX2UB"] prefix = "CPLEX2UB_" has_cplex_opts = True elif key.startswith("CPLEXDUAL_"): section = config_lines["CPLEXDUAL"] prefix = "CPLEXDUAL_" has_cplex_opts = True elif key.startswith("CPLEX2DUAL_"): section = config_lines["CPLEX2DUAL"] prefix = "CPLEX2DUAL_" has_cplex_opts = True else: section = config_lines[None] prefix = "" val = self.options[key] line = "%s %s\n" % (key.replace(prefix,'',1), val) section.append(line) with open(config_filename, "r+") as f: f.seek(append_pos) for line in config_lines[None]: f.write(line) if has_cplex_opts: f.write("\nCPLEXBEGIN\n") for line in config_lines["CPLEX"]: f.write(line) for section_key in ('CPLEXEEV', 'CPLEXLB', 'CPLEX2LB', 'CPLEXUB', 'CPLEX2UB', 'CPLEXDUAL', 'CPLEX2DUAL'): section = config_lines[section_key] if len(section) > 0: f.write(section_key+"\n") for line in section: f.write(line) f.write("CPLEXEND\n") f.write("\n\nEND\n") def _read_solution(self, sp, symbols_filename, info_filename, solution_filename): """Parses a DDSIP solution file.""" # parse the symbol map symbol_map = {} with open(symbols_filename) as f: for line in f: lp_symbol, scenario_tree_id = line.strip().split() symbol_map[lp_symbol] = scenario_tree_id # # Xhat # try: xhat = {} with open(solution_filename, 'r') as f: line = f.readline() while line.strip() != "1. Best Solution": line = f.readline() line = f.readline() assert line.startswith("Variable name Value") line = f.readline() assert line.startswith("-----------------------------------") line = f.readline().strip() while line != "": line = line.split() varlabel, varsol = line xhat[symbol_map[varlabel]] = float(varsol) line = f.readline().strip() except (IOError, OSError): logger.warn( "Exception encountered while parsing ddsip " "solution file '%s':\n%s'" % (solution_filename, traceback.format_exc())) xhat = None # # Objective, bound, status, etc. # results = SPSolverResults() results.solver.status_code = None results.status = None results.solver.status = None results.solver.termination_condition = None results.solver.message = None results.solver.time = None try: with open(info_filename, 'r') as f: line = f.readline() while True: if line.startswith("Total CPU time:"): break line = f.readline() if line == '': # Unexpected file format or the solve failed logger.warn( "Unexpected ddsip info file format. No " "status information will be returned") return xhat, results line = f.readline().strip() while (line == "") or \ (line == "----------------------------------------------------------------------------------------"): line = f.readline().strip() if line.startswith("EEV"): results.eev = float(line.split()[1]) line = f.readline().strip() if line.startswith("VSS"): results.vss = float(line.split()[1]) line = f.readline().strip() assert line.startswith("EVPI") line = line.split() results.evpi = float(line[1]) line = f.readline().strip() assert line == "" line = f.readline().strip() assert line == "----------------------------------------------------------------------------------------" line = f.readline().strip() line = line.split() assert len(line) == 4 assert line[0] == "Status" results.solver.status_code = int(line[1]) assert line[2] == "Time" results.solver.time = float(line[3]) (results.status, results.solver.status, results.solver.termination_condition, results.solver.message) = \ _ddsip_status_map.get(results.solver.status_code, (SolutionStatus.unknown, SolverStatus.unknown, TerminationCondition.unknown, None)) line = f.readline().strip() line = line.split() assert len(line) == 6 assert line[0] == "Upper" assert line[3] == "Tree" results.tree_depth = int(line[5]) line = f.readline().strip() line = line.split() assert len(line) == 2 assert line[0] == "Nodes" results.nodes = int(line[1]) line = f.readline().strip() while line != "----------------------------------------------------------------------------------------": line = f.readline().strip() line = f.readline().strip() assert line.startswith("Best Value") results.objective = float(line.split()[2]) # NOTE: I think DDSIP refers to the "bound" as # "Lower Bound", even when the objective # is being maximized. line = f.readline().strip() assert line.startswith("Lower Bound") results.bound = float(line.split()[2]) except (IOError, OSError): logger.warn( "Exception encountered while parsing ddsip " "info file '%s':\n%s'" % (info_filename, traceback.format_exc())) return xhat, results def runddsip_register_options(options=None): if options is None: options = PySPConfigBlock() DDSIPSolver.register_options(options) ScenarioTreeManagerFactory.register_options(options) safe_register_common_option(options, "verbose") safe_register_common_option(options, "disable_gc") safe_register_common_option(options, "profile") safe_register_common_option(options, "traceback") safe_register_common_option(options, "output_scenario_tree_solution") safe_register_common_option(options, "keep_solver_files") safe_register_common_option(options, "output_solver_log") safe_register_common_option(options, "symbolic_solver_labels") # used to populate the implicit DDSIP options safe_register_unique_option( options, "solver_options", PySPConfigValue( (), domain=_domain_tuple_of_str_or_dict, description=( "Unregistered solver options that will be passed " "to DDSIP via the config file (e.g., NODELIM=4, " "CPLEX_1067=1). This option can be used multiple " "times from the command line to specify more " "than one DDSIP option." ), doc=None, visibility=0), ap_kwds={'action': 'append'}, ap_group=_ddsip_group_label) return options def runddsip(options): """ Construct a senario tree manager and solve it with the DDSIP solver. """ start_time = time.time() with ScenarioTreeManagerFactory(options) as sp: sp.initialize() print("") print("Running DDSIP solver for stochastic " "programming problems") ddsip = DDSIPSolver() # add the implicit ddsip options solver_options = options.solver_options if len(solver_options) > 0: if type(solver_options) is tuple: for name_val in solver_options: assert "=" in name_val name, val = name_val.split("=") ddsip.options[name.strip()] = val.strip() else: for key, val in solver_options.items(): ddsip.options[key] = val ddsip_options = ddsip.extract_user_options_to_dict(options, sparse=True) results = ddsip.solve( sp, options=ddsip_options, output_solver_log=options.output_solver_log, keep_solver_files=options.keep_solver_files, symbolic_solver_labels=options.symbolic_solver_labels) xhat = results.xhat del results.xhat print("") print(results) if options.output_scenario_tree_solution: print("") sp.scenario_tree.snapshotSolutionFromScenarios() sp.scenario_tree.pprintSolution() sp.scenario_tree.pprintCosts() print("") print("Total execution time=%.2f seconds" % (time.time() - start_time)) return 0 # # the main driver routine # def main(args=None): # # Top-level command that executes everything # # # Import plugins # import pyomo.environ # # Parse command-line options. # try: options = parse_command_line( args, runddsip_register_options, prog='runddsip', description=( """Optimize a stochastic program using the DDSIP solver.""" )) except SystemExit as _exc: # the parser throws a system exit if "-h" is specified # - catch it to exit gracefully. return _exc.code return launch_command(runddsip, options, error_label="runddsip: ", disable_gc=options.disable_gc, profile_count=options.profile, traceback=options.traceback) SPSolverFactory.register_solver("ddsip", DDSIPSolver) if __name__ == "__main__": sys.exit(main())