# _________________________________________________________________________ # # 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 logging import traceback import pyutilib.subprocess 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 ScenarioTreeManagerFactory import pyomo.pysp.convert.ddsip 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())