Changeset 2068 for coopr.pysp/stable


Ignore:
Timestamp:
Dec 30, 2009 12:43:59 AM (11 years ago)
Author:
wehart
Message:

Merged revisions 1952-2067 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pysp/trunk

........

r1956 | jwatson | 2009-12-02 17:56:53 -0700 (Wed, 02 Dec 2009) | 3 lines


Added --scenario-solver-options and --ef-solver-options options to the "runph" script.

........

r1957 | dlwoodr | 2009-12-03 14:17:35 -0700 (Thu, 03 Dec 2009) | 2 lines


Documentation updates for pysp

........

r1974 | wehart | 2009-12-06 17:20:56 -0700 (Sun, 06 Dec 2009) | 2 lines


Updating PyPI categories

........

r1978 | jwatson | 2009-12-10 21:29:33 -0700 (Thu, 10 Dec 2009) | 3 lines


Eliminated exception-handling logic when loading user-defined extension modules in PH. The range of exceptions is too large, and for debugging purposes, it is more useful to see the raw trace output.

........

r1979 | jwatson | 2009-12-10 22:23:17 -0700 (Thu, 10 Dec 2009) | 5 lines


Biggest enhancement: The efwriter command-line script now has the option to output a CVaR-weighted objective term. Not extensively tested, but behaves sane on a number of small test cases.


Improved exception handling and error diagnostics in both the runph and efwriter scripts.

........

r1985 | jwatson | 2009-12-12 10:45:17 -0700 (Sat, 12 Dec 2009) | 3 lines


Modified PH to only use warm-starts if a solver has the capability!

........

r1998 | jwatson | 2009-12-13 15:17:58 -0700 (Sun, 13 Dec 2009) | 3 lines


Changed references to _component to active_component.

........

r2026 | wehart | 2009-12-21 23:27:06 -0700 (Mon, 21 Dec 2009) | 2 lines


Attempting to update PH. I'm not sure if this works, since I don't know how to test PH.

........

r2029 | jwatson | 2009-12-22 09:52:21 -0700 (Tue, 22 Dec 2009) | 3 lines


Some fixes to the ef writer based on Bill's recent changes to _initialize_constraint.

........

r2035 | jwatson | 2009-12-22 21:10:32 -0700 (Tue, 22 Dec 2009) | 3 lines


Added --scenario-mipgap option to PH script. Added _mipgap attribute to PH object. This attribute is mirrored to the solver plugin at the initiation of each iteration, after any PH extensions have the opportunity to provide a new value to the attribute. It is currently made use of by the WW PH extension.

........

r2037 | dlwoodr | 2009-12-23 14:38:43 -0700 (Wed, 23 Dec 2009) | 2 lines


add this example from Fernando

........

r2038 | dlwoodr | 2009-12-23 14:46:56 -0700 (Wed, 23 Dec 2009) | 3 lines


finish the job: we can now duplicate Fernando's example

........

r2039 | jwatson | 2009-12-23 15:13:54 -0700 (Wed, 23 Dec 2009) | 3 lines


Missed fix with new constraint initialization syntax in PH linearization.

........

r2058 | jwatson | 2009-12-29 10:57:58 -0700 (Tue, 29 Dec 2009) | 3 lines


Missed some _initialize_constraint function calls within the PySP EF writer during the recent switch to the corresponding "add" method.

........

r2059 | jwatson | 2009-12-29 10:58:34 -0700 (Tue, 29 Dec 2009) | 3 lines


Enabling garbage collection by default in PH.

........

r2060 | jwatson | 2009-12-29 10:59:05 -0700 (Tue, 29 Dec 2009) | 3 lines


Elimnating some debug output.

........

r2061 | jwatson | 2009-12-29 11:07:47 -0700 (Tue, 29 Dec 2009) | 3 lines


Fixing some option documentation in PH.

........

r2062 | jwatson | 2009-12-29 12:00:37 -0700 (Tue, 29 Dec 2009) | 3 lines


Added ef-mipgap option to PH scripts.

........

Location:
coopr.pysp/stable/2.1
Files:
18 edited
26 copied

Legend:

Unmodified
Added
Removed
  • coopr.pysp/stable/2.1

  • coopr.pysp/stable/2.1/coopr/pysp/asynchph.py

    r1768 r2068  
    280280      self._converger = converger
    281281
    282       model_objective = model._component[Objective]
     282      model_objective = model.active_components(Objective)
    283283      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
    284284
     
    532532         for scenario in self._scenario_tree._scenarios:
    533533            instance = self._instances[scenario._name]
    534             for objective_name in instance._component[Objective]:
    535                objective = instance._component[Objective][objective_name]
     534            for objective_name in instance.active_components(Objective):
     535               objective = instance.active_components(Objective)objective_name]
    536536               # TBD: I don't know how to deal with objective arrays, so I'll assume they aren't there
    537537               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
     
    692692      for instance_name, instance in self._instances.items():
    693693         
    694          objective_name = instance._component[Objective].keys()[0]         
    695          objective = instance._component[Objective][objective_name]         
     694         objective_name = instance.active_components(Objective).keys()[0]         
     695         objective = instance.active_components(Objective)[objective_name]         
    696696         objective_expression = objective._data[None].expr # TBD: we don't deal with indexed expressions (not even sure how to interpret them)
    697697         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
     
    711711
    712712                  w_parameter_name = "PHWEIGHT_"+variable_name
    713                   w_parameter = instance._component[param._ParamBase][w_parameter_name]
     713                  w_parameter = instance.active_components(Param)[w_parameter_name]
    714714                 
    715715                  average_parameter_name = "PHAVG_"+variable_name
    716                   average_parameter = instance._component[param._ParamBase][average_parameter_name]
     716                  average_parameter = instance.active_components(Param)[average_parameter_name]
    717717
    718718                  rho_parameter_name = "PHRHO_"+variable_name
    719                   rho_parameter = instance._component[param._ParamBase][rho_parameter_name]
     719                  rho_parameter = instance.active_components(Param)[rho_parameter_name]
    720720
    721721                  for index in variable_indices:
    722722
    723                      instance_variable = instance._component[var._VarBase][variable_name][index]
     723                     instance_variable = instance.active_components(Var)[variable_name][index]
    724724                     
    725725                     if (instance_variable.status is not VarStatus.unused) and (instance_variable.fixed is False):
     
    733733         # we'll be covered.
    734734         objective_expression.simplify(instance)
    735          instance._component[Objective][objective_name]._data[None].expr = objective_expression
    736          instance._component[Objective][objective_name]._quad_subexpr = quad_expression
     735         instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
     736         instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
    737737
    738738   def iteration_k_plus_solves(self):
     
    803803
    804804        if self._verbose == True:
    805            for objective_name in instance._component[Objective]:
    806               objective = instance._component[Objective][objective_name]
     805           for objective_name in instance.active_components(Objective):
     806              objective = instance.active_components(Objective)[objective_name]
    807807              print "%20s       %18.4f     %14.4f" % (scenario_name, ph_objective_value, 0.0)
    808808
  • coopr.pysp/stable/2.1/coopr/pysp/ef.py

    r1774 r2068  
    2121
    2222   # check the master model first.
    23    for var in master_model._component[_VarBase].values():
     23   for var in master_model.active_components(Var).values():
    2424      if isinstance(var.domain, BooleanSet):
    2525         return True
     
    2828   for scenario_name in scenario_instances.keys():
    2929      scenario_instance = scenario_instances[scenario_name]
    30       for var in scenario_instance._component[_VarBase].values():
     30      for var in scenario_instance.active_components(Var).values():
    3131         if isinstance(var.domain, BooleanSet):
    3232            return True
     
    4040
    4141   # check the master model first.
    42    for var in master_model._component[_VarBase].values():
     42   for var in master_model.active_components(Var).values():
    4343      if isinstance(var.domain, IntegerSet):
    4444         return True
     
    4747   for scenario_name in scenario_instances.keys():
    4848      scenario_instance = scenario_instances[scenario_name]
    49       for var in scenario_instance._component[_VarBase].values():
     49      for var in scenario_instance.active_components(Var).values():
    5050         if isinstance(var.domain, IntegerSet):
    5151            return True
     
    5454
    5555# the main extensive-form writer routine - including read of scenarios/etc.
    56 def write_ef_from_scratch(model_directory, instance_directory, output_filename):
     56def write_ef_from_scratch(model_directory, instance_directory, output_filename, \
     57                          generate_weighted_cvar, cvar_weight, risk_alpha):
    5758
    5859   start_time = time.time()
     
    6667   #### INITIALIZATION ############################################################################
    6768   ################################################################################################
     69
     70   #
     71   # validate cvar options, if specified.
     72   #
     73   if generate_weighted_cvar is True:
     74      if cvar_weight <= 0.0:
     75         raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)
     76      if (risk_alpha <= 0.0) or (risk_alpha >= 1.0):
     77         raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)
     78
     79      print "Writing CVaR weighted objective"
     80      print "CVaR term weight="+str(cvar_weight)
     81      print "CVaR alpha="+str(risk_alpha)
     82      print ""
    6883   
    6984   #
     
    261276                        new_constraint = Constraint(name=new_constraint_name)
    262277                        new_expr = master_variable[index] - scenario_variable[index]
    263                         new_constraint._initialize_constraint(None, (0.0, new_expr, 0.0), master_binding_instance, new_constraint_name)
     278                        new_constraint.add(None, (0.0, new_expr, 0.0))
    264279                        new_constraint._model = master_binding_instance
    265280                        setattr(master_binding_instance, new_constraint_name, new_constraint)
     
    272287            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
    273288
     289   # if we're generating the weighted CVaR objective term, create the corresponding variable and
     290   # the master CVaR eta variable.
     291   if generate_weighted_cvar is True:
     292      root_node = scenario_tree._stages[0]._tree_nodes[0]
     293     
     294      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
     295      cvar_cost_variable = Var(name=cvar_cost_variable_name)
     296      setattr(master_binding_instance, cvar_cost_variable_name, cvar_cost_variable)
     297      cvar_cost_variable.construct()
     298
     299      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     300      cvar_eta_variable = Var(name=cvar_eta_variable_name)
     301      setattr(master_binding_instance, cvar_eta_variable_name, cvar_eta_variable)     
     302      cvar_eta_variable.construct()
     303
    274304   master_binding_instance.presolve()
    275305
     
    308338            new_constraint = Constraint(name=new_constraint_name)
    309339            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
    310             new_constraint._initialize_constraint(None, (0.0, new_expr, 0.0), master_binding_instance, new_constraint_name)
     340            new_constraint.add(None, (0.0, new_expr, 0.0))
    311341            new_constraint._model = master_binding_instance
    312342            setattr(master_binding_instance, new_constraint_name, new_constraint)
     
    338368         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
    339369         new_constraint = Constraint(name=new_constraint_name)
    340          new_constraint._initialize_constraint(None, (0.0, constraint_expr, 0.0), master_binding_instance, new_constraint_name)
     370         new_constraint.add(None, (0.0, constraint_expr, 0.0))
    341371         new_constraint._model = master_binding_instance                     
    342372         setattr(master_binding_instance, new_constraint_name, new_constraint)
     
    345375
    346376            an_instance = instances[instances.keys()[0]]
    347             an_objective = an_instance._component[Objective]
     377            an_objective = an_instance.active_components(Objective)
    348378            opt_sense = an_objective[an_objective.keys()[0]].sense
    349379
     380            opt_expression = node_expected_cost_variable
     381
     382            if generate_weighted_cvar is True:
     383               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
     384               cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
     385               opt_expression += cvar_weight * cvar_cost_variable
     386
    350387            new_objective = Objective(name="MASTER", sense=opt_sense)
    351             new_objective._data[None].expr = node_expected_cost_variable
     388            new_objective._data[None].expr = opt_expression
    352389            setattr(master_binding_instance, "MASTER", new_objective)
     390
     391   # CVaR requires the addition of a variable per scenario to represent the cost excess,
     392   # and a constraint to compute the cost excess relative to eta. we also replicate (following
     393   # what we do for node cost variables) an eta variable for each scenario instance, and
     394   # require equality with the master eta variable via constraints.
     395   if generate_weighted_cvar is True:
     396     
     397      root_node = scenario_tree._stages[0]._tree_nodes[0]
     398
     399      master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     400      master_cvar_eta_variable = getattr(master_binding_instance, master_cvar_eta_variable_name)
     401     
     402      for scenario_name in instances.keys():
     403         scenario_instance = instances[scenario_name]
     404
     405         # unique names are required because the presolve isn't
     406         # aware of the "owning" models for variables.
     407         cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
     408         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
     409         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
     410         cvar_excess_variable.construct()
     411
     412         cvar_eta_variable_name = "CVAR_ETA"
     413         cvar_eta_variable = Var(name=cvar_eta_variable_name)
     414         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
     415         cvar_eta_variable.construct()
     416
     417         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
     418         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
     419         compute_excess_expression = cvar_excess_variable
     420         for node in scenario_tree._scenario_map[scenario_name]._node_list:
     421            (cost_variable, cost_variable_idx) = node._stage._cost_variable
     422            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
     423         compute_excess_expression += cvar_eta_variable
     424         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
     425         compute_excess_constraint._model = scenario_instance
     426         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
     427
     428         eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
     429         eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
     430         eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
     431         eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
     432         eta_equality_constraint._model = master_binding_instance
     433         setattr(master_binding_instance, eta_equality_constraint_name, eta_equality_constraint)
     434
     435      # add the constraint to compute the master CVaR variable value. iterate
     436      # over scenario instances to create the expected excess component first.
     437      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
     438      cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
     439      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     440      cvar_eta_variable = getattr(master_binding_instance, cvar_eta_variable_name)
     441     
     442      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
     443     
     444      for scenario_name in instances.keys():
     445         scenario_instance = instances[scenario_name]
     446         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
     447
     448         scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
     449         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
     450
     451         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
     452
     453      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
     454      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
     455      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
     456      compute_cvar_cost_constraint._model = master_binding_instance
     457      setattr(master_binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
     458
     459   # after mucking with instances, presolve to collect terms required prior to output.         
     460   # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
     461   for scenario_name in instances.keys():
     462      scenario_instance = instances[scenario_name]   
     463      scenario_instance.presolve()         
    353464
    354465   master_binding_instance.presolve()
     
    572683                        new_constraint = Constraint(name=new_constraint_name)
    573684                        new_expr = master_variable[index] - scenario_variable[index]
    574                         new_constraint._initialize_constraint(None, (0.0, new_expr, 0.0), master_binding_instance, new_constraint_name)
     685                        new_constraint.add(None, (0.0, new_expr, 0.0))
    575686                        new_constraint._model = master_binding_instance
    576687                        setattr(master_binding_instance, new_constraint_name, new_constraint)
     
    619730            new_constraint = Constraint(name=new_constraint_name)
    620731            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
    621             new_constraint._initialize_constraint(None, (0.0, new_expr, 0.0), master_binding_instance, new_constraint_name)
     732            new_constraint.add(None, (0.0, new_expr, 0.0))
    622733            new_constraint._model = master_binding_instance
    623734            setattr(master_binding_instance, new_constraint_name, new_constraint)
     
    649760         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
    650761         new_constraint = Constraint(name=new_constraint_name)
    651          new_constraint._initialize_constraint(None, (0.0, constraint_expr, 0.0), master_binding_instance, new_constraint_name)
     762         new_constraint.add(None, (0.0, constraint_expr, 0.0))
    652763         new_constraint._model = master_binding_instance                     
    653764         setattr(master_binding_instance, new_constraint_name, new_constraint)
     
    656767
    657768            an_instance = instances[instances.keys()[0]]
    658             an_objective = an_instance._component[Objective]
     769            an_objective = an_instance.active_components(Objective)
    659770            opt_sense = an_objective[an_objective.keys()[0]].sense
    660771
  • coopr.pysp/stable/2.1/coopr/pysp/ef_writer_script.py

    r1768 r2068  
    4040                  type="string",
    4141                  default=".")
     42parser.add_option("--generate-weighted-cvar",
     43                  help="Add a weighted CVaR term to the primary objective",
     44                  action="store_true",
     45                  dest="generate_weighted_cvar",
     46                  default=False)
     47parser.add_option("--cvar-weight",
     48                  help="The weight associated with the CVaR term in the risk-weighted objective formulation.",
     49                  action="store",
     50                  dest="cvar_weight",
     51                  type="float",
     52                  default=1.0)
     53parser.add_option("--risk-alpha",
     54                  help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.",
     55                  action="store",
     56                  dest="risk_alpha",
     57                  type="float",
     58                  default=0.95)
    4259parser.add_option("--output-file",
    4360                  help="Specify the name of the extensive form output file",
     
    5572def run_ef_writer(options, args):
    5673
    57    write_ef_from_scratch(options.model_directory, options.instance_directory, options.output_file)
     74   # if the user enabled the addition of the weighted cvar term to the objective,
     75   # then validate the associated parameters.
     76   generate_weighted_cvar = False
     77   cvar_weight = None
     78   risk_alpha = None
     79
     80   if options.generate_weighted_cvar is True:
     81
     82      generate_weighted_cvar = True
     83      cvar_weight = options.cvar_weight
     84      risk_alpha = options.risk_alpha
     85
     86   write_ef_from_scratch(options.model_directory, options.instance_directory, options.output_file, \
     87                         generate_weighted_cvar, cvar_weight, risk_alpha)
    5888
    5989   return
  • coopr.pysp/stable/2.1/coopr/pysp/ph.py

    r1952 r2068  
    475475               # this is a quadratic penalty term - the lower bound is 0!
    476476               new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
    477                new_penalty_term_variable._construct(instance,None)                                             
     477               new_penalty_term_variable.construct()
    478478               setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
    479479               self._instance_augmented_attributes[instance_name].append(new_penalty_term_variable_name)
     
    594594      self._integer_tolerance = 0.00001
    595595
     596      # PH maintains a mipgap that is applied to each scenario solve that is performed.
     597      # this attribute can be changed by PH extensions, and the change will be applied
     598      # on all subsequent solves - until it is modified again. the default is None,
     599      # indicating unassigned.
     600      self._mipgap = None
     601
     602      # we only store these temporarily...
     603      scenario_solver_options = None
     604
    596605      # process the keyword options
    597606      for key in kwds.keys():
     
    607616            self._solver_type = kwds[key]
    608617         elif key == "solver_manager":
    609             self._solver_manager_type = kwds[key]           
     618            self._solver_manager_type = kwds[key]
     619         elif key == "scenario_solver_options":
     620            scenario_solver_options = kwds[key]
     621         elif key == "scenario_mipgap":
     622            self._mipgap = kwds[key]           
    610623         elif key == "keep_solver_files":
    611624            self._keep_solver_files = kwds[key]
     
    637650         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
    638651      if self._rho <= 0.0:
    639          raise ValueError, "Value of rho paramter in PH must be non-zero positive; value specified=" + `self._rho`
     652         raise ValueError, "Value of the rho parameter in PH must be non-zero positive; value specified=" + `self._rho`
     653      if (self._mipgap is not None) and ((self._mipgap < 0.0) or (self._mipgap > 1.0)):
     654         raise ValueError, "Value of the mipgap parameter in PH must be on the unit interval; value specified=" + `self._mipgap`
    640655
    641656      # validate the linearization (number of pieces) and breakpoint distribution parameters.
     
    669684      if self._keep_solver_files is True:
    670685         self._solver.keepFiles = True
     686      if len(scenario_solver_options) > 0:
     687         if self._verbose is True:
     688            print "Initializing scenario sub-problem solver with options="+str(scenario_solver_options)
     689         self._solver.set_options("".join(scenario_solver_options))
    671690
    672691      # construct the solver manager.
     
    731750      self._converger = converger
    732751
    733       model_objective = model._component[Objective]
     752      model_objective = model.active_components(Objective)
    734753      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
    735754
     
    844863      self._original_objective_expression = {}
    845864      for instance_name, instance in self._instances.items():
    846          objective_name = instance._component[Objective].keys()[0]
    847          self._original_objective_expression[instance_name] = instance._component[Objective][objective_name]._data[None].expr
     865         objective_name = instance.active_components(Objective).keys()[0]
     866         self._original_objective_expression[instance_name] = instance.active_components(Objective)[objective_name]._data[None].expr
    848867
    849868      # cache the number of discrete and continuous variables in the master instance. this value
     
    889908
    890909      solve_start_time = time.time()
     910
     911      # STEP 0: set up all global solver options.
     912      self._solver.mipgap = self._mipgap
    891913
    892914      # STEP 1: queue up the solves for all scenario sub-problems and
     
    958980         for scenario in self._scenario_tree._scenarios:
    959981            instance = self._instances[scenario._name]
    960             for objective_name in instance._component[Objective]:
    961                objective = instance._component[Objective][objective_name]
     982            for objective_name in instance.active_components(Objective):
     983               objective = instance.active_components(Objective)[objective_name]
    962984               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
    963985         print "------------------------------------------------"
     
    10801102      for instance_name, instance in self._instances.items():
    10811103         
    1082          objective_name = instance._component[Objective].keys()[0]         
    1083          objective = instance._component[Objective][objective_name]
     1104         objective_name = instance.active_components(Objective).keys()[0]         
     1105         objective = instance.active_components(Objective)[objective_name]
    10841106         # clone the objective, because we're going to augment (repeatedly) the original objective.
    10851107         objective_expression = self._original_objective_expression[instance_name].clone()
     
    11021124
    11031125               w_parameter_name = "PHWEIGHT_"+variable_name
    1104                w_parameter = instance._component[param._ParamBase][w_parameter_name]
     1126               w_parameter = instance.active_components(Param)[w_parameter_name]
    11051127                 
    11061128               average_parameter_name = "PHAVG_"+variable_name
    1107                average_parameter = instance._component[param._ParamBase][average_parameter_name]
     1129               average_parameter = instance.active_components(Param)[average_parameter_name]
    11081130
    11091131               rho_parameter_name = "PHRHO_"+variable_name
    1110                rho_parameter = instance._component[param._ParamBase][rho_parameter_name]
     1132               rho_parameter = instance.active_components(Param)[rho_parameter_name]
    11111133
    11121134               blend_parameter_name = "PHBLEND_"+variable_name
    1113                blend_parameter = instance._component[param._ParamBase][blend_parameter_name]
     1135               blend_parameter = instance.active_components(Param)[blend_parameter_name]
    11141136
    11151137               node_min_parameter = variable_tree_node._minimums[variable_name]
     
    11191141               if self._linearize_nonbinary_penalty_terms > 0:
    11201142                  quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
    1121                   quad_penalty_term_variable = instance._component[var._VarBase][quad_penalty_term_variable_name]
    1122 
    1123                instance_variable = instance._component[var._VarBase][variable_name]
     1143                  quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name]
     1144
     1145               instance_variable = instance.active_components(Var)[variable_name]
    11241146
    11251147               for index in variable_indices:
     
    11831205                                 ub = x.ub()
    11841206
     1207                              if x.lb == x.ub:
     1208                                 print "***WARNING - LB EQUALS UB"
     1209                                 print "VALUE=",x.lb
     1210                                 print "VARIABLE="+variable_name+indexToString(index)
     1211
    11851212                              node_min = node_min_parameter[index]()
    11861213                              node_max = node_max_parameter[index]()
     
    12111238                                 piece_constraint.model = instance
    12121239                                 piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index])
    1213                                  piece_constraint._initialize_constraint(None, piece_expression, instance, piece_constraint_name)
     1240                                 piece_constraint.add(None, piece_expression)
    12141241                                 setattr(instance, piece_constraint_name, piece_constraint)                                   
    12151242
     
    12221249         # we'll be covered.
    12231250         objective_expression.simplify(instance)
    1224          instance._component[Objective][objective_name]._data[None].expr = objective_expression
     1251         instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
    12251252         # if we are linearizing everything, then nothing will appear in the quadratic expression -
    12261253         # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
    12271254         # be properly generated.
    12281255         if quad_expression != 0.0:
    1229            instance._component[Objective][objective_name]._quad_subexpr = quad_expression
     1256           instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
    12301257
    12311258   def iteration_k_solve(self):
     
    12391266
    12401267     solve_start_time = time.time()
     1268
     1269     # STEP 0: set up all global solver options.
     1270     self._solver.mipgap = self._mipgap     
    12411271
    12421272     # STEP 1: queue up the solves for all scenario sub-problems and
     
    12611291        # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
    12621292        new_action_handle = None
    1263         if self._disable_warmstarts is False:
     1293        if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
    12641294           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
    12651295        else:
     
    13141344        for scenario in self._scenario_tree._scenarios:
    13151345           instance = self._instances[scenario._name]
    1316            for objective_name in instance._component[Objective]:
    1317               objective = instance._component[Objective][objective_name]
     1346           for objective_name in instance.active_components(Objective):
     1347              objective = instance.active_components(Objective)[objective_name]
    13181348              print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
    13191349
  • coopr.pysp/stable/2.1/coopr/pysp/ph_script.py

    r1952 r2068  
    7070                  type="string",
    7171                  default="serial")
     72parser.add_option("--scenario-solver-options",
     73                  help="Solver options for all PH scenario sub-problems",
     74                  action="append",
     75                  dest="scenario_solver_options",
     76                  type="string",
     77                  default=[])
     78parser.add_option("--scenario-mipgap",
     79                  help="Specifies the mipgap for all PH scenario sub-problems",
     80                  action="store",
     81                  dest="scenario_mipgap",
     82                  type="float",
     83                  default=None)
     84parser.add_option("--ef-solver-options",
     85                  help="Solver options for the extension form problem",
     86                  action="append",
     87                  dest="ef_solver_options",
     88                  type="string",
     89                  default=[])
     90parser.add_option("--ef-mipgap",
     91                  help="Specifies the mipgap for the EF solve",
     92                  action="store",
     93                  dest="ef_mipgap",
     94                  type="float",
     95                  default=None)
    7296parser.add_option("--max-iterations",
    7397                  help="The maximal number of PH iterations. Default is 100.",
     
    223247                  default=0)
    224248parser.add_option("--restore-from-checkpoint",
    225                   help="The name of the checkpoint file from which PH should be initialized. Default is "", indicating no checkpoint restoration",
     249                  help="The name of the checkpoint file from which PH should be initialized. Default is \"\", indicating no checkpoint restoration",
    226250                  action="store",
    227251                  dest="restore_from_checkpoint",
     
    235259                  default=0)
    236260parser.add_option("--enable-gc",
    237                   help="Enable the python garbage collecter.",
     261                  help="Enable the python garbage collecter. Default is True.",
    238262                  action="store_true",
    239263                  dest="enable_gc",
    240                   default=False)
     264                  default=True)
    241265
    242266
     
    424448
    425449      if options.user_defined_extension is not None:
    426          print "User-defined PH extension specified - module name="+options.user_defined_extension
    427          try:
    428             __import__(options.user_defined_extension)
    429          except:
    430             raise RuntimeError, "Failed to load user-defined module="+options.user_defined_extension
     450         print "Trying to import user-defined PH extension module="+options.user_defined_extension
     451         # JPW removed the exception handling logic, as the module importer
     452         # can raise a broad array of exceptions.
     453         __import__(options.user_defined_extension)
     454         print "Module successfully loaded"
    431455
    432456      #
     
    452476                              solver=options.solver_type, \
    453477                              solver_manager=options.solver_manager_type, \
     478                              scenario_solver_options=options.scenario_solver_options, \
     479                              scenario_mipgap=options.scenario_mipgap, \
    454480                              keep_solver_files=options.keep_solver_files, \
    455481                              output_solver_log=options.output_solver_logs, \
     
    513539      if ef_solver is None:
    514540         raise ValueError, "Failed to create solver of type="+options.solver_type+" for use in extensive form solve"
     541      if len(options.ef_solver_options) > 0:
     542         print "Initializing ef solver with options="+str(options.ef_solver_options)         
     543         ef_solver.set_options("".join(options.ef_solver_options))
     544      if options.ef_mipgap is not None:
     545         if (options.ef_mipgap < 0.0) or (options.ef_mipgap > 1.0):
     546            raise ValueError, "Value of the mipgap parameter for the EF solve must be on the unit interval; value specified=" + `options.ef_mipgap`
     547         else:
     548            ef_solver.mipgap = options.ef_mipgap
    515549
    516550      ef_solver_manager = SolverManagerFactory(options.solver_manager_type)
     
    523557      ef_results = ef_solver_manager.wait_for(ef_action_handle)
    524558      print "Extensive form solve results:"
    525       print ef_results
    526      
     559      ef_results.write(num=1)
    527560
    528561def run(args=None):
  • coopr.pysp/stable/2.1/coopr/pysp/scenariotree.py

    r1768 r2068  
    5151      # IMPT: This implicitly assumes convergence across the scenarios - if not, garbage results.
    5252      instance = scenario_instance_map[self._scenarios[0]._name]
    53       my_cost = instance._component[var._VarBase][self._stage._cost_variable[0].name][self._stage._cost_variable[1]]()
     53      my_cost = instance.active_components(Var)[self._stage._cost_variable[0].name][self._stage._cost_variable[1]]()
    5454      child_cost = 0.0
    5555      for child in self._children:
     
    294294
    295295               # validate that the variable exists and extract the reference.
    296                if variable_name not in self._reference_instance._component[var._VarBase]:
     296               if variable_name not in self._reference_instance.active_components(Var):
    297297                  raise ValueError, "Variable=" + variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    298                variable = self._reference_instance._component[var._VarBase][variable_name]               
     298               variable = self._reference_instance.active_components(Var)[variable_name]               
    299299
    300300               # extract all "real", i.e., fully specified, indices matching the index template.
     
    311311
    312312               # verify that the variable exists.
    313                if variable_string not in self._reference_instance._component[var._VarBase].keys():
     313               if variable_string not in self._reference_instance.active_components(Var).keys():
    314314                  raise RuntimeError, "Unknown variable=" + variable_string + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    315315
    316                variable = self._reference_instance._component[var._VarBase][variable_string]
     316               variable = self._reference_instance.active_components(Var)[variable_string]
    317317
    318318               # 9/14/2009 - now forcing the user to explicit specify the full
     
    349349
    350350            # validate that the variable exists and extract the reference.
    351             if cost_variable_name not in self._reference_instance._component[var._VarBase]:
     351            if cost_variable_name not in self._reference_instance.active_components(Var):
    352352               raise ValueError, "Variable=" + cost_variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    353             cost_variable = self._reference_instance._component[var._VarBase][cost_variable_name]               
     353            cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]               
    354354
    355355            # extract all "real", i.e., fully specified, indices matching the index template.
     
    367367
    368368            # validate that the variable exists and extract the reference
    369             if cost_variable_name not in self._reference_instance._component[var._VarBase]:
     369            if cost_variable_name not in self._reference_instance.active_components(Var):
    370370               raise ValueError, "Cost variable=" + cost_variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    371             cost_variable = self._reference_instance._component[var._VarBase][cost_variable_name]
     371            cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]
    372372           
    373373         # store the validated info.
     
    582582         aggregate_cost = 0.0
    583583         for stage in self._stages:
    584             instance_cost_variable = instance._component[var._VarBase][stage._cost_variable[0].name][stage._cost_variable[1]]()
     584            instance_cost_variable = instance.active_components(Var)[stage._cost_variable[0].name][stage._cost_variable[1]]()
    585585            print "\tStage=%20s     Cost=%10.4f" % (stage._name, instance_cost_variable)
    586586            aggregate_cost += instance_cost_variable
  • coopr.pysp/stable/2.1/coopr/pysp/wwphextension.py

    r1952 r2068  
    11#  _________________________________________________________________________
    22#
    3 coopr: a common optimization python repository
    4 copyright (c) 2009 sandia corporation.
     3Coopr: A COmmon Optimization Python Repository
     4Copyright (c) 2009 Sandia Corporation.
    55#  this software is distributed under the bsd license.
    66#  under the terms of contract de-ac04-94al85000 with sandia corporation,
     
    100100               
    101101               # verify that the root variable exists and grab it.
    102                if variable_name not in reference_instance._component[var._VarBase].keys():
     102               if variable_name not in reference_instance.active_components(Var).keys():
    103103                  raise RuntimeError, "Unknown variable="+variable_name+" referenced in ww ph extension suffix file="+self._suffix_filename
    104                variable = reference_instance._component[var._VarBase][variable_name]
     104               variable = reference_instance.active_components(Var)[variable_name]
    105105
    106106               # extract all "real", i.e., fully specified, indices matching the index template.
     
    132132
    133133               # verify that the variable exists.
    134                if variable_string not in reference_instance._component[var._VarBase].keys():
     134               if variable_string not in reference_instance.active_components(Var).keys():
    135135                  raise RuntimeError, "Unknown variable="+variable_string+" referenced in ww ph extension suffix file="+self._suffix_filename
    136136
    137                variable = reference_instance._component[var._VarBase][variable_string]
     137               variable = reference_instance.active_components(Var)[variable_string]
    138138
    139139               # 9/14/2009 - now forcing the user to explicit specify the full
     
    323323      if self.Iteration0MipGap > 0.0:
    324324         print "Setting mipgap to "+str(self.Iteration0MipGap)
    325          ph._solver.options.mip = "tolerances mipgap " + str(self.Iteration0MipGap)
     325         ph._mipgap = self.Iteration0MipGap
    326326         
    327327#==================================================
     
    371371                              # update convergence prior to checking for fixing.
    372372                              self._int_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
    373                               attrvariable = ph._model_instance._component[var._VarBase][variable_name][index]
     373                              attrvariable = ph._model_instance.active_components(Var)[variable_name][index]
    374374                              if hasattr(attrvariable, 'Iter0FixIfConvergedAtLB'):
    375375                                 lb = getattr(attrvariable, 'Iter0FixIfConvergedAtLB')
     
    409409         gap = self.InitialMipGap
    410410         print "Setting mipgap to "+str(gap)
    411          ph._solver.options.mip = "tolerances mipgap " + str(gap)
     411         ph._mipgap = gap
    412412
    413413#==================================================
     
    458458
    459459                              # now check on permissions to converge to various placed (e.g., lb is lb permission)
    460                               attrvariable = ph._model_instance._component[var._VarBase][variable_name][index]
     460                              attrvariable = ph._model_instance.active_components(Var)[variable_name][index]
    461461                              if hasattr(attrvariable, 'FixWhenItersConvergedAtLB'):
    462462                                 lb = getattr(attrvariable, 'FixWhenItersConvergedAtLB')
     
    531531         else:
    532532            print "Setting mipgap to "+str(gap)
    533          ph._solver.options.mip = "tolerances mipgap " + str(gap)
     533         ph._mipgap = gap
    534534
    535535
     
    748748
    749749         # verify that the root variable exists and grab it.
    750          if variable_name not in reference_instance._component[var._VarBase].keys():
     750         if variable_name not in reference_instance.active_components(Var).keys():
    751751            raise RuntimeError, "Unknown variable="+variable_name+" referenced while slamming. "
    752          variable = reference_instance._component[var._VarBase][variable_name]
     752         variable = reference_instance.active_components(Var)[variable_name]
    753753
    754754         didone = False;   # did we find at least one node to slam in?
  • coopr.pysp/stable/2.1/doc/pysp/pyspbody.tex

    r1830 r2068  
    1717
    1818The pysp package extends the pyomo modeling language to support multi-stage stochastic programs with
    19 enumerated scenarios. Pyomo and pysp are Python version 2.6 programs. The underlying algorithm in pysp
    20 is based on Progressive Hedging (PH) \cite{RockafellarWets}, which progressively computes {\em weights} corresponding to each variable to force convergence. In order to
     19enumerated scenarios. Pyomo and pysp are Python version 2.6 programs.  In order to
    2120specify a program, the user must provide a reference model and a scenario tree.
    2221
     
    3231\end{verbatim}
    3332but notice that there are two dashes before the word ``help.'' Command line arguments are summarized in Section~\ref{cmdargsec}.
     33
     34The underlying algorithm in pysp
     35is based on Progressive Hedging (PH) \cite{RockafellarWets}, which decomposes the problem into sub-problems, one for
     36each scenario.
     37The algorithm progressively computes {\em weights} corresponding to each variable to force convergence
     38and also makes use of a {\em proximal} term that provides a penalty for the squared deviation from the
     39mean solution from the last PH iteration.
    3440
    3541\subsection{Reference Model}
     
    140146Note that both PH extensions (WW PH and your own) can co-exist; however, the WW plugin will be invoked first.
    141147 
     148\item \verb|--scenario-solver-options| The options are specified just as in pyomo, e.g., \verb|--scenario-solver-options="mip_tolerances_mipgap=0.2"| to set the mipgap for all scenario sub-problem solves to 20\% for the CPLEX solver. The options are specified in a quote deliminted string that is passed to the sub-problem solver.
     149 Whatever options specified are persistent across all solves. 
     150\item \verb|--ef-solver-options| The options are specified just as in pyomo, e.g., \verb|--scenario-solver-options="mip_tolerances_mipgap=0.2"| to set the mipgap for all scenario sub-problem solves to 20\% for the CPLEX solver. The options are specified in a quote deliminted string that is passed to the EF problem solver.
    142151  \item \verb|--write-ef|\\            Upon termination, write the extensive form of the model - accounting for all fixed variables.
    143152  \item \verb|--solve-ef|\\            Following write of the extensive form model, solve it.
     
    146155  \item \verb|--suppress-continuous-variable-output|\\
    147156                        Eliminate PH-related output involving continuous variables. Default: no output.
    148   \item --keep-solver-files\\   Retain temporary input and output files for scenario sub-problem solves.  Default: files not kept.
    149   \item --output-solver-logs\\  Output solver logs during scenario sub-problem solves. Default: no output.
    150   \item --output-ef-solver-log\\
     157  \item \verb|--keep-solver-files|\\   Retain temporary input and output files for scenario sub-problem solves.  Default: files not kept.
     158  \item \verb|--output-solver-logs|\\  Output solver logs during scenario sub-problem solves. Default: no output.
     159  \item \verb|--output-ef-solver-log|\\
    151160                        Output solver log during the extensive form solve. Default: no output.
    152   \item --output-solver-results\\
     161  \item \verb|--output-solver-results|\\
    153162                        Output solutions obtained after each scenario sub-problem solve. Default: no output.
    154   \item --output-times\\        Output timing statistics for various PH components. Default: no output.
     163  \item \verb|--output-times|\\        Output timing statistics for various PH components. Default: no output.
    155164  \item \verb|--disable-warmstarts|\\
    156165                  Disable warm-start of scenario sub-problem solves in PH iterations >= 1.
     
    162171                  Do not linearize PH objective terms involving binary decision variables.
    163172                  Default=False (i.e., the proximal term for binary variables is linearized by default; this can have some impact on the relaxations during the branch and bound solution process).
    164   \item \verb|: --linearize-nonbinary-penalty-terms|\\
    165                   Approximate the PH quadratic term for continuous variables with a piece-wise linear function.
    166                   Default=False.
     173  \item \verb|--linearize-nonbinary-penalty-terms|=BPTS\\
     174                  Approximate the PH quadratic term for non-binary variables with a piece-wise linear function. The argument
     175                  BPTS gives the number of breakpoints in the linear approximation. The default=0. Reasonable non-zero
     176                  values are usually in the range of 3 to 7. Note that if a breakpoint would be very close to a variable
     177                  bound, then the break point is ommited. IMPORTANT: this option requires that all variables have bounds that
     178                  are either established in the reference model or by code specfied using the bounds-cfgfile command line option. See Section~\ref{LinearSec}
     179                  for more information about linearizing the proximal term.
     180  \item \verb|--breakpoint-strategy|=BREAKPOINT\_STRATEGY
     181                        Specify the strategy to distribute breakpoints on the
     182                        [lb, ub] interval of each variable when linearizing. 0
     183                        indicates uniform distribution. 1 indicates
     184                        breakpoints at the node min and max, uniformly in-
     185                        between. 2 indicates more aggressive concentration of
     186                        breakpoints near the observed node min/max.
     187  \item \verb|--bounds-cfgfile|=BOUNDS\_CFGFILE\\
     188                  The argument BOUNDS\_CFGFILE specifies the name of an executable pyomo file that sets bounds. The devault is that there is no file.
     189                  When specified, the code in this file is executed after the initialization of scenario data so the
     190                  bounds can be based on data from all scenarios. The config subdirectory of the farmer example contains a simple example of
     191                  such a file (boundsetter.cfg).
    167192  \item \verb|--checkpoint-interval|\\
    168                   The number of iterations between writing of a checkpoint file. Default is 0, indicating never..
     193                  The number of iterations between writing of a checkpoint file. Default is 0, indicating never.
    169194  \item \verb|--restore-from-checkpoint|\\
    170195                  The name of the checkpoint file from which PH should be initialized. Default is not to restore from a checkpoint.
     
    344369
    345370\begin{verbatim}
    346 pyomo\examples\pysp\farmer
    347 \end{verbatim}
    348 where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
    349 
    350 The data is in a subdirectory called SCENARIODATA.
    351 
    352 To
     371coopr\examples\pysp\farmer
     372\end{verbatim}
     373
     374The data is in a subdirectory called scenariodata and the model is in the models subdirectory. To
    353375invoke PH for this problem, connect to this farmer directory and use the command:
    354376
     
    379401
    380402\begin{verbatim}
    381 pyomo\packages\coopr\examples\pysp\sizes
    382 \end{verbatim}
    383 where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
     403coopr\packages\coopr\examples\pysp\sizes
     404\end{verbatim}
    384405
    385406The data for a three scenario version is in a subdirectory called SIZES3 and a ten scenario
     
    415436
    416437\begin{verbatim}
    417 pyomo\packages\coopr\examples\pysp\forestry
    418 \end{verbatim}
    419 where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
     438coopr\packages\coopr\examples\pysp\forestry
     439\end{verbatim}
     440
    420441
    421442There are two families of instances: ``Chile'' and ``Davis,'' each with four stages and eighteen scenarios. This is also a small
     
    10101031To run this example, connect to the sizes directory, which is something like:
    10111032\begin{verbatim}
    1012 pyomo\packages\coopr\examples\pysp\sizes
    1013 \end{verbatim}
    1014 where the directory \verb|pyomo| might be \verb|pyomodist| or something else depending on how you installed pyomo and pysp.
     1033coopr\examples\pysp\sizes
     1034\end{verbatim}
    10151035
    10161036Then use
     
    10251045time. It is provided to illustrate the features of the extensions, not to illustrate why you might need them. Much larger instances are are required for that.
    10261046
     1047\section{Linearizing the Proximal Term \label{LinearSec}}
     1048
     1049\subsection{Introduction}
     1050
     1051For a decision variable $x$ the proximal term added to the objective function for each scenario subproblem at PH iteration $k$ is
     1052$$
     1053\left(x - \bar{x}^{(k-1)}\right)^{2}
     1054$$
     1055where $\bar{x}^{(k-1)}$ is the average over the scenarios from the last iteration. Expanding the square reveals that the only
     1056quadratic term is $x^{2}$. For binary variables, this is equal to $x$, although this is not the case when a relaxation is solved.
     1057For binary variables, the default behaviour is to replace the quadratic term with $x$, but an option allows the quadratic to be
     1058retained because this can effect the subproblem solution time due to the use of the quadratic term when the relaxations are solved
     1059as part of the branch and bound process.
     1060
     1061For non-binary variables an option exists to replace the $x^2$ term with a piece-wise approximation and the number of breakpoints in the
     1062approximation is under user control. This can have a significant effect on CPU time required to solve sub-problems because the presence of
     1063the quadratic term increases the solution times significantly. However, linearization results in a time-quality tradeoff because increasing
     1064the number of breakpoints increases the fidelity but each breakpoint introduces another (unseen) binary variable so solution times are generally
     1065increased.
     1066
     1067A few strategies for placing the breakpoints are supported as command a line options:
     1068\verb|--breakpoint-strategy|=BREAKPOINT\_STRATEGY. A value of 0
     1069indicates uniform distribution. 1 indicates
     1070breakpoints at the min and max for the node in the scenario tree, uniformly in-
     1071between. A value of 2 indicates more aggressive concentration of
     1072breakpoints near the observed node min/max.
     1073
     1074Upper and lower bounds for variables must be specified when the linearization option is chosen. This can be done by specifying bounds in the
     1075reference model or by using the bounds-cfgfile command line option. It is, of course, best to use meaningful bounds provided
     1076in the reference model; however, the modeller must be careful not use estimated bounds that are too tight since that could preclude
     1077an optimal (or even a feasible) solution to the overall stochastic problem even though it might not cut off any optimal solutions for the
     1078particular scenario. The use of a bounds cfgfile is an advanced topic, but enables the modeller to use bounds that cannot create
     1079infeasibilities.
     1080
     1081\subsection{Bounds-cfgfile}
     1082
     1083Using a bounds cfgfile is an advanced topic. The modeler is writing python/pyomo code that will be executed by the
     1084ph.py file that is the core of the PH algorithm. The first executable line in a bounds file is typically:
     1085
     1086\begin{verbatim}
     1087model_instance = self._model_instance
     1088\end{verbatim}
     1089
     1090This establishes a local variable called ``model\_instance'' that can be used to refer to the model (of course, a different name can be
     1091used, such as MODINS). The object ``self'' on the right hand side of this assignment refers to the core PH object. The model, in turn contains the parameters and
     1092variables defined in the ReferenceModel.py file that can be accessed by name. For example, with the Farmer model, the cfg file
     1093sets the uppter bound on DevotedAcreage to be value of the paramter TOTAL\_ACREAGE at intialization (since this is Python, the
     1094parentheses after TOTAL\_ACREAGE cause the value in TOTAL\_ACREAGE to be assigned rather than the name of the parameter object:
     1095
     1096\begin{verbatim}
     1097# only need to set upper bounds on first-stage variables, i.e., those being blended.
     1098
     1099model_instance = self._model_instance
     1100
     1101# the values supplied to the method
     1102upper_bound = float(model_instance.TOTAL_ACREAGE())
     1103
     1104for index in model_instance.CROPS:
     1105   # the lower and upper bounds are expected to be floats, and trigger an exception if not.
     1106   self.setVariableBoundsAllScenarios("DevotedAcreage", index, 0.0, upper_bound)
     1107\end{verbatim}
     1108
     1109The same thing could be accomplished by setting the upper bound in the model file, but it does serve as a simple example of a bounds cfg file.
     1110
     1111%\subsection{Sizes Example}
     1112
     1113
    10271114\section{Solving sub-problems in Parallel \label{parallelsec}}
    10281115
     
    10331120Here are the commands in order:
    10341121\begin{enumerate}
    1035 \item On the master: \verb|pyro-ns|
     1122\item On the master: \verb|coopr-ns|
    10361123\item On the master: \verb|dispatch_srvr|
    10371124\item On each slave: \verb|pyro_mip_server|
    10381125\item On the master: \verb|runph ... --solver-manager=pyro ...|
    10391126\end{enumerate}
    1040 Note that the command \verb|pryo-ns| and the argument \verb|solver-manger| have a dash in the middle, while
     1127Note that the command \verb|coopr-ns| and the argument \verb|solver-manger| have a dash in the middle, while
    10411128the commands \verb|dispatch_srvr| and \verb|pyro_mip_server| have underscores. The first three commands launch processes
    10421129that have no internal mechanism for termination; i.e., they will be terminated only if they crash or if they are
    1043 killed by an external process. It is common to launch these processes with output redirection, such as \verb|pyro-ns >& pyrons.log|. The
     1130killed by an external process. It is common to launch these processes with output redirection, such as \verb|coopr-ns >& cooprns.log|. The
    10441131\verb|runph| command is a normal runph command with the usual arguments with the additional specification that subproblem solves should
    10451132be directed to the pyro solver manager.
  • coopr.pysp/stable/2.1/examples/pysp/forestry/config/boundsetter.cfg

    r1952 r2068  
    88
    99# formula taken from Fernando's AMPL model. Not sure why the time period is restricted to the first year.
    10 umax = sum([model_instance.A[h] * model_instance.a[h,"Ano1"] for h in model_instance.HarvestCells])
    11 self.setVariableBoundsAllIndicesAllScenarios("f", 0.0, umax())
     10#umax = sum([model_instance.A[h] * model_instance.a[h,"Ano1"] for h in model_instance.HarvestCells])
     11#self.setVariableBoundsAllIndicesAllScenarios("f", 0.0, umax())
    1212
    1313# can be eliminated once Pyomo supports parameters in variable bounds declarations.
     
    1919         instance.z[e, t].setlb(lb)
    2020         instance.z[e, t].setub(ub)
     21
     22for instance_name, instance in self._instances.items():
     23   for t in instance.Times:
     24      lb = 0
     25      ub = instance.Zub[t]()
     26      for (i,j) in instance.AllRoads:
     27         instance.f[i,j,t].setub(ub)
     28         instance.f[i,j,t].setlb(lb)
  • coopr.pysp/stable/2.1/examples/pysp/forestry/config/rhosetter.cfg

    r1754 r2068  
    44model_instance = self._model_instance
    55
    6 MyRhoFactor = 0.1
     6MyRhoFactorDelta = 0.001
     7MyRhoFactorGamma = 0.0001
     8MyRhoFactorF = 0.1
     9MyRhoFactorZ = 0.01
    710
    8 for t in model_instance.Times:
     11for t in model_instance.Times[:-1]:
    912
    1013   for h in model_instance.HarvestCells:
    11       self.setRhoAllScenarios( model_instance.delta[h,t], model_instance.P[h,t] * model_instance.A[h] * MyRhoFactor)
     14      self.setRhoAllScenarios( model_instance.delta[h,t] , model_instance.P[h,t] * model_instance.A[h] * MyRhoFactorDelta + MyRhoFactorDelta * model_instance.a[h, t] * model_instance.yr[t] * model_instance.A[h] * sum([ model_instance.Q[o, t] for o in model_instance.COriginNodeForCell[h]]) )
    1215      # .a[h, t] * .A[h] * sum([ model_instance.Q[o, t] for o in model.OriginNodes for h in model.HCellsForOrigin[o]])
     16      # model_instance.a[h, t] * model_instance.A[h] * sum([ model_instance.Q[o, t] for o in model_instance.COriginNodeForCell[h]])
    1317
    1418   for (i,j) in model_instance.PotentialRoads:
    15       self.setRhoAllScenarios( model_instance.gamma[i,j,t], model_instance.C[i,j,t] * MyRhoFactor)
     19      self.setRhoAllScenarios( model_instance.gamma[i,j,t], model_instance.C[i,j,t] * MyRhoFactorGamma)
    1620
    1721   for (i,j) in model_instance.AllRoads:
    18       self.setRhoAllScenarios( model_instance.f[i,j,t], model_instance.D[i,j,t] * MyRhoFactor)
     22      self.setRhoAllScenarios( model_instance.f[i,j,t], model_instance.D[i,j,t] * MyRhoFactorF)
    1923
    2024   for e in model_instance.ExitNodes:
    2125   #for e in model.ExitNodes:
    22       self.setRhoAllScenarios( model_instance.z[e,t],  model_instance.R[e,t] * MyRhoFactor)
     26      self.setRhoAllScenarios( model_instance.z[e,t],  model_instance.R[e,t] * MyRhoFactorZ)
    2327
  • coopr.pysp/stable/2.1/examples/pysp/forestry/config/wwph-nb.suffixes

    r1859 r2068  
    22
    33# The default should be zero at iterations 0, apart from Gamma
    4 gamma[*,*,Ano1] Iter0FixIfConvergedAtUB 1
    5 gamma[*,*,Ano2] Iter0FixIfConvergedAtUB 1
    6 gamma[*,*,Ano3] Iter0FixIfConvergedAtUB 0
    7 gamma[*,*,Ano4] Iter0FixIfConvergedAtUB 0
     4#gamma[*,*,Ano1] Iter0FixIfConvergedAtUB 0
     5#gamma[*,*,Ano2] Iter0FixIfConvergedAtUB 0
     6#gamma[*,*,Ano3] Iter0FixIfConvergedAtUB 0
     7#gamma[*,*,Ano4] Iter0FixIfConvergedAtUB 0
    88
    9 gamma[*,*,*] Iter0FixIfConvergedAtUB 1
     9gamma[*,*,*] Iter0FixIfConvergedAtUB 0
     10gamma[*,*,*] Iter0FixIfConvergedAtLB 0
     11
     12delta[*,*] Iter0FixIfConvergedAtUB 0
     13delta[*,*] Iter0FixIfConvergedAtLB 0
    1014
    1115# if you want all indices, you can just use the variable name (no [*,*]).
    12 gamma[*,*,*] FixWhenItersConvergedAtLB 8
    13 gamma[*,*,*] FixWhenItersConvergedAtUB 4
     16gamma[*,*,Ano1] FixWhenItersConvergedAtLB 15
     17gamma[*,*,Ano1] FixWhenItersConvergedAtUB 10
     18gamma[*,*,Ano2] FixWhenItersConvergedAtLB 16
     19gamma[*,*,Ano2] FixWhenItersConvergedAtUB 11
     20gamma[*,*,Ano3] FixWhenItersConvergedAtLB 17
     21gamma[*,*,Ano3] FixWhenItersConvergedAtUB 12
    1422
    15 delta[*,*] FixWhenItersConvergedAtLB 10
    16 delta[*,*] FixWhenItersConvergedAtUB 10
     23delta[*,Ano1] FixWhenItersConvergedAtLB 15
     24delta[*,Ano1] FixWhenItersConvergedAtUB 10
     25delta[*,Ano2] FixWhenItersConvergedAtLB 16
     26delta[*,Ano2] FixWhenItersConvergedAtUB 11
     27delta[*,Ano3] FixWhenItersConvergedAtLB 17
     28delta[*,Ano3] FixWhenItersConvergedAtUB 12
    1729
    1830gamma[*,*,*] CanSlamToLB False
    1931gamma[*,*,*] CanSlamToMin False
    2032gamma[*,*,*] CanSlamToAnywhere False
    21 gamma[*,*,*] CanSlamToMax True
    22 gamma[*,*,*] CanSlamToUB True
     33gamma[*,*,*] CanSlamToMax False
     34gamma[*,*,*] CanSlamToUB False
    2335
    2436delta[*,*] CanSlamToLB False
    2537delta[*,*] CanSlamToMin False
    26 delta[*,*] CanSlamToAnywhere True
     38delta[*,*] CanSlamToAnywhere False
    2739delta[*,*] CanSlamToMax False
    2840delta[*,*] CanSlamToUB False
    2941
    3042# higher priority means slam these first
    31 gamma[*,*,Ano1] SlammingPriority 40
    32 gamma[*,*,Ano2] SlammingPriority 30
    33 gamma[*,*,Ano3] SlammingPriority 20
    34 gamma[*,*,Ano4] SlammingPriority 10
    35      
    36 delta[*,*] SlammingPriority 2
    37 
    38 
     43#gamma[*,*,Ano1] SlammingPriority 40
     44#gamma[*,*,Ano2] SlammingPriority 30
     45#gamma[*,*,Ano3] SlammingPriority 20
     46#gamma[*,*,Ano4] SlammingPriority 10
     47#delta[*,*] SlammingPriority 2
  • coopr.pysp/stable/2.1/examples/pysp/forestry/config/wwph.cfg

    r1754 r2068  
    11self.Iter0FixIfConvergedAtLB = 0
    2 self.Iter0FixIfConvergedAtUB = 1
     2self.Iter0FixIfConvergedAtUB = 0
    33self.Iter0FixIfConvergedAtNB = 0
    44
    5 self.FixWhenItersConvergedAtLB = 10 
    6 self.FixWhenItersConvergedAtUB = 10
    7 self.FixWhenItersConvergedAtNB = 10
     5#self.FixWhenItersConvergedAtLB = 10 
     6#self.FixWhenItersConvergedAtUB = 10
     7#self.FixWhenItersConvergedAtNB = 10
    88
    99self.Iteration0MipGap = 0.025
    10 self.InitialMipGap = 0.25
    11 self.FinalMipGap = 0.0005  # will be above this
     10self.InitialMipGap = 0.05
     11self.FinalMipGap = 0.005  # will be above this
    1212
    1313self.fix_continuous_variables = False
    1414self.StartBlendingContIter = 2
     15self.DisableCycleDetection = True
  • coopr.pysp/stable/2.1/examples/pysp/forestry/models-nb-yr/ReferenceModel.py

    r1759 r2068  
    4646model.HCellsForOrigin = Set(model.OriginNodes, within=model.HarvestCells)
    4747
     48# HERE'S MY BUG!!!!!!!!!!!!!!!!!!!
     49# origin nodes for cell
     50#def origin_node_for_harvest_cell_rule(h, model):
     51#   for o in model.OriginNodes:
     52#      if h in model.HCellsForOrigin[o]:
     53#         print "h,o=",h,o
     54#         return o
     55#         #return o()
     56#         #return str(o)
     57#model.COriginNodeForCellA = Set(model.HarvestCells, within=model.OriginNodes, rule=origin_node_for_harvest_cell_rule)
     58# SO I PRECALCULATED IT
     59model.COriginNodeForCell = Set(model.HarvestCells, within=model.OriginNodes)
     60
    4861# existing roads
    4962model.ExistingRoads = Set(within=model.Nodes*model.Nodes)
     
    6679model.Aminus = Set(model.Nodes, within=model.AllRoads, rule=aminus_arc_set_rule)
    6780
     81# derived set for the "no isolated roads can be built" constraint
     82# if you would build the road, then one connecting road should be built around you. And because the roads its isolated, you need to sum over potential roads around you
     83def connected_potential_roads_rule(k, l, model):
     84   return  [ (i, j) for (i, j) in model.PotentialRoads if (i==k and j!=l) or (i==l and j!=k) or (j==k and i!=l) or (j==l and i!=k) ]
     85model.ConnectedPotentialRoads = Set( model.Nodes, model.Nodes, within=model.PotentialRoads, rule=connected_potential_roads_rule)
     86# a PotentialRoad is not isolated if it connects to an ExitNode or an Existing Road
     87def potential_roads_connected_to_existing_roads_rule(model):
     88   set=[]
     89   for (k,l) in model.PotentialRoads:
     90      #if not (k in model.ExitNodes and l in model.OriginNodes) and not (l in model.ExitNodes and k in model.OriginNodes):
     91      if not (k in model.ExitNodes) and not (l in model.ExitNodes):
     92         for x in model.Nodes:
     93            if x!=k and x!=l:
     94               #print "k,l,x=",k,l,x
     95               if (x,k) in model.AllRoads or (x,l) in model.AllRoads or (k,x) in model.AllRoads or (l,x) in model.AllRoads:
     96                  #if not (k in model.OriginNodes and (l,x) in model.ExistingRoads) and not (l in model.OriginNodes and (k,x) in model.ExistingRoads):
     97                  if ((x,k) in model.ExistingRoads) or ((x,l) in model.ExistingRoads) or ((k,x) in model.ExistingRoads) or ((l,x) in model.ExistingRoads):
     98                     if (k,l) not in set:
     99                        set.append((k,l))
     100                        #print "set=",set
     101   return set
     102model.PotentialRoadsConnectedToExistingRoads = Set(within=model.PotentialRoads, rule=potential_roads_connected_to_existing_roads_rule)
     103# isolated PotentialRoads
     104model.DisconnectedPotentialRoads = model.PotentialRoads - model.PotentialRoadsConnectedToExistingRoads
     105
     106# derived set for the "no isolated lots can be harvest" constraint
     107# if a isolated lot is harvested then one potential road accesing it must be built around
     108def potential_roads_accesing_lot_h(h, model):
     109   return  [ (i, j) for (i, j) in model.PotentialRoads if (i in model.COriginNodeForCell[h]) or (j in model.COriginNodeForCell[h]) ]
     110model.PotentialRoadsAccesingLotH = Set( model.HarvestCells, within=model.PotentialRoads, rule=potential_roads_accesing_lot_h )
     111# a isolated lot is not connected to an ExistingRoad
     112def lots_connected_to_existing_roads(model):
     113   set=[]
     114   for h in model.HarvestCells:
     115      for (i,j) in model.ExistingRoads:
     116         if (i in model.COriginNodeForCell[h]) or (j in model.COriginNodeForCell[h]):
     117            if h not in set:
     118               set.append(h)
     119   return set
     120model.LotsConnectedToExistingRoads = Set(within=model.HarvestCells, rule=lots_connected_to_existing_roads )
     121# isolated lots or HarvestCells
     122model.LotsNotConnectedToExistingRoads =  model.HarvestCells - model.LotsConnectedToExistingRoads
     123
     124
    68125#
    69126# Deterministic Parameters
     
    78135# flow capacity (smallest big M)
    79136def Umax_init(model):
    80    return (sum([1.1 * model.a[h, t] * model.A[h] for h in model.HarvestCells for t in model.Times])/4)()
     137   return (sum([model.a[h, t] * model.A[h] for h in model.HarvestCells for t in model.Times])/2)()
    81138model.Umax = Param(initialize=Umax_init)
    82139
     
    104161model.Zub = Param(model.Times)
    105162
    106 # yield ratio of the forest
     163# Yield Ratio
    107164model.yr = Param(model.Times)
    108165
     
    192249
    193250model.EnforceOneHarvest = Constraint(model.HarvestCells, rule=one_harvest)
     251
     252# no isolated roads can be built
     253def road_to_road_trigger(t, k, l, model):
     254   return ( None, sum([ model.gamma[ k, l, tt] for tt in model.Times if tt <= t ]) - sum([ model.gamma[ i, j, ttt] for ttt in model.Times if ttt <= t for (i,j) in model.ConnectedPotentialRoads[k,l] ]), 0.0 )
     255
     256model.EnforceRoadToRoadTrigger = Constraint(model.Times, model.DisconnectedPotentialRoads, rule=road_to_road_trigger)
     257
     258# no isolated lots can be harvest
     259def lot_to_road_trigger(t, h, model):
     260   return ( None, sum([ model.delta[ h, tt] for tt in model.Times if tt <= t ]) - sum([ model.gamma[ i, j, ttt] for ttt in model.Times if ttt <= t for (i,j) in model.PotentialRoadsAccesingLotH[h] ]), 0.0 )
     261
     262model.EnforceLotToRoadTrigger = Constraint(model.Times, model.LotsNotConnectedToExistingRoads, rule=lot_to_road_trigger)
    194263
    195264#
  • coopr.pysp/stable/2.1/scripts/efwriter

    r1806 r2068  
    5151    coopr.pysp.ef_writer_script.run()
    5252except ValueError, str:
     53    print "VALUE ERROR:"   
    5354    print str
    5455except IOError, str:
     56    print "IO ERROR:"   
    5557    print str
     58except RuntimeError, str:
     59    print "RUN-TIME ERROR:"   
     60    print str   
    5661except pyutilib.common.ApplicationError, str:
     62    print "APPLICATION ERROR:"   
    5763    print str
    5864except:
  • coopr.pysp/stable/2.1/scripts/runph

    r1806 r2068  
    5959    print "APPLICATION ERROR:"
    6060    print str
     61except RuntimeError, str:
     62    print "RUN-TIME ERROR:"   
     63    print str       
    6164except:
    6265   print "Encountered unhandled exception"
  • coopr.pysp/stable/2.1/setup.py

    r1952 r2068  
    6161            'Programming Language :: Unix Shell',
    6262            'Topic :: Scientific/Engineering :: Mathematics',
    63             'Topic :: Software Development :: Libraries :: Python Modules'
     63            'Topic :: Software Development :: Libraries :: Python Modules',
     64            'Topic :: Scientific/Engineering'
    6465        ],
    6566      packages=packages,
Note: See TracChangeset for help on using the changeset viewer.