diff --git a/doc/src/extensions.rst b/doc/src/extensions.rst index 5274ee7e9..cc223f5d5 100644 --- a/doc/src/extensions.rst +++ b/doc/src/extensions.rst @@ -116,6 +116,23 @@ constructor or in the hub dictionary under ``opt_kwargs`` as the There is an example of the function in the sizes example (``_rho_setter``). +SepRho +^^^^^^ + +Set per variable rho values using the "SEP" algorithm from + +Progressive hedging innovations for a class of stochastic mixed-integer resource allocation problems +Jean-Paul Watson, David L. Woodruff, Compu Management Science, 2011 +DOI 10.1007/s10287-010-0125-4 + +One can additional specify a multiplier on the computed value (default = 1.0). +If the cost coefficient on a non-anticipative variable is 0, the default rho value is used instead. + +CoeffRho +^^^^^^^^ + +Set per variable rho values proportional to the cost coefficient on each non-anticipative variable, +with an optional multiplier (default = 1.0). If the coefficient is 0, the default rho value is used instead. wtracker_extension ^^^^^^^^^^^^^^^^^^ diff --git a/examples/afew.py b/examples/afew.py index 5bb5520ae..89fb52a7b 100644 --- a/examples/afew.py +++ b/examples/afew.py @@ -47,7 +47,7 @@ def do_one(dirname, progname, np, argstring): # for farmer, the first arg is num_scens and is required do_one("farmer", "farmer_cylinders.py", 3, "--num-scens=3 --bundles-per-rank=0 --max-iterations=50 " - "--default-rho=1 --display-convergence-detail " + "--default-rho=1 --sep-rho --display-convergence-detail " "--solver-name={} --xhatshuffle --lagrangian --use-norm-rho-updater".format(solver_name)) do_one("farmer", "farmer_lshapedhub.py", 2, "--num-scens=3 --bundles-per-rank=0 --max-iterations=50 " @@ -58,7 +58,7 @@ def do_one(dirname, progname, np, argstring): 4, "--num-scens=3 --bundles-per-rank=0 --max-iterations=5 " "--iter0-mipgap=0.01 --iterk-mipgap=0.001 --linearize-proximal-terms " - " --xhatshuffle --lagrangian --fwph " + " --xhatshuffle --lagrangian --fwph --smoothing " "--default-rho=1 --solver-name={} --display-progress".format(solver_name)) do_one("hydro", "hydro_cylinders_pysp.py", 3, diff --git a/examples/farmer/farmer_cylinders.py b/examples/farmer/farmer_cylinders.py index 11d516676..acec3fcb2 100644 --- a/examples/farmer/farmer_cylinders.py +++ b/examples/farmer/farmer_cylinders.py @@ -45,6 +45,8 @@ def _parse_args(): cfg.wxbar_read_write_args() cfg.tracking_args() cfg.reduced_costs_args() + cfg.sep_rho_args() + cfg.coeff_rho_args() cfg.add_to_config("crops_mult", description="There will be 3x this many crops (default 1)", domain=int, @@ -125,6 +127,11 @@ def main(): if cfg.reduced_costs: vanilla.add_reduced_costs_fixer(hub_dict, cfg) + if cfg.sep_rho: + vanilla.add_sep_rho(hub_dict, cfg) + if cfg.coeff_rho: + vanilla.add_coeff_rho(hub_dict, cfg) + # FWPH spoke if cfg.fwph: fw_spoke = vanilla.fwph_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs) @@ -145,6 +152,10 @@ def main(): ph_ob_spoke = vanilla.ph_ob_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs, rho_setter = rho_setter) + if cfg.sep_rho: + vanilla.add_sep_rho(ph_ob_spoke, cfg) + if cfg.coeff_rho: + vanilla.add_coeff_rho(ph_ob_spoke, cfg) # xhat looper bound spoke if cfg.xhatlooper: diff --git a/examples/run_all.py b/examples/run_all.py index e36d4851d..88014f238 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -286,11 +286,10 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): 4, "--instance-name=sslp_15_45_10 --bundles-per-rank=2 " "--max-iterations=5 --default-rho=1 " - "--subgradient --xhatshuffle --fwph " + "--subgradient --xhatshuffle --fwph --coeff-rho " "--linearize-proximal-terms " "--rel-gap=0.0 " "--solver-name={} --fwph-stop-check-tol 0.01".format(solver_name)) - do_one("sslp", "sslp_cylinders.py", 3, diff --git a/examples/sslp/sslp_cylinders.py b/examples/sslp/sslp_cylinders.py index 638004240..3397e6b5e 100644 --- a/examples/sslp/sslp_cylinders.py +++ b/examples/sslp/sslp_cylinders.py @@ -32,6 +32,7 @@ def _parse_args(): cfg.xhatshuffle_args() cfg.subgradient_args() cfg.reduced_costs_args() + cfg.coeff_rho_args() cfg.parse_command_line("sslp_cylinders") return cfg @@ -87,9 +88,15 @@ def main(): if reduced_costs: vanilla.add_reduced_costs_fixer(hub_dict, cfg) + if cfg.coeff_rho: + vanilla.add_coeff_rho(hub_dict, cfg) + # FWPH spoke if fwph: fw_spoke = vanilla.fwph_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs) + # Need to fix FWPH to support extensions + # if cfg.coeff_rho: + # vanilla.add_coeff_rho(fw_spoke, cfg) # Standard Lagrangian bound spoke if lagrangian: @@ -100,6 +107,8 @@ def main(): subgradient_spoke = vanilla.subgradient_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs, rho_setter = None) + if cfg.coeff_rho: + vanilla.add_coeff_rho(subgradient_spoke, cfg) # xhat looper bound spoke if xhatlooper: diff --git a/mpisppy/extensions/coeff_rho.py b/mpisppy/extensions/coeff_rho.py new file mode 100644 index 000000000..3aa687923 --- /dev/null +++ b/mpisppy/extensions/coeff_rho.py @@ -0,0 +1,40 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### + +import mpisppy.extensions.extension + +from mpisppy.utils.sputils import nonant_cost_coeffs + + +class CoeffRho(mpisppy.extensions.extension.Extension): + """ + Determine rho as a linear function of the objective coefficient + """ + + def __init__(self, ph): + self.ph = ph + self.multiplier = 1.0 + if ( + "coeff_rho_options" in ph.options + and "multiplier" in ph.options["coeff_rho_options"] + ): + self.multiplier = ph.options["coeff_rho_options"]["multiplier"] + + def post_iter0(self): + for s in self.ph.local_scenarios.values(): + cc = nonant_cost_coeffs(s) + for ndn_i, rho in s._mpisppy_model.rho.items(): + if cc[ndn_i] != 0: + rho._value = abs(cc[ndn_i]) * self.multiplier + # if self.ph.cylinder_rank==0: + # nv = s._mpisppy_data.nonant_indices[ndn_i] # var_data object + # print(ndn_i,nv.getname(),cc[ndn_i],rho._value) + + if self.ph.cylinder_rank == 0: + print("Rho values updated by CoeffRho Extension") diff --git a/mpisppy/extensions/reduced_costs_fixer.py b/mpisppy/extensions/reduced_costs_fixer.py index 6b36f053a..be4625624 100644 --- a/mpisppy/extensions/reduced_costs_fixer.py +++ b/mpisppy/extensions/reduced_costs_fixer.py @@ -196,7 +196,7 @@ def reduced_costs_fixing(self, reduced_costs): if self.opt.cylinder_rank == 0 and self.verbose: print("All reduced costs are nan, heuristic fixing will not be applied") return - + # compute the quantile target abs_reduced_costs = np.abs(reduced_costs) diff --git a/mpisppy/extensions/sep_rho.py b/mpisppy/extensions/sep_rho.py new file mode 100644 index 000000000..2fea016ff --- /dev/null +++ b/mpisppy/extensions/sep_rho.py @@ -0,0 +1,170 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### + +import mpisppy.extensions.extension +import numpy as np +import mpisppy.MPI as MPI + +from mpisppy.utils.sputils import nonant_cost_coeffs + + +class SepRho(mpisppy.extensions.extension.Extension): + """ + Rho determination algorithm "SEP" from + Progressive hedging innovations for a class of stochastic mixed-integer + resource allocation problems + Jean-Paul Watson, David L. Woodruff, Compu Management Science, 2011 + DOI 10.1007/s10287-010-0125-4 + """ + + def __init__(self, ph): + self.ph = ph + + self.multiplier = 1.0 + + if ( + "sep_rho_options" in ph.options + and "multiplier" in ph.options["sep_rho_options"] + ): + self.multiplier = ph.options["sep_rho_options"]["multiplier"] + + def _compute_primal_residual_norm(self, ph): + local_nodenames = [] + local_primal_residuals = {} + global_primal_residuals = {} + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + for node in s._mpisppy_node_list: + if node.name not in local_nodenames: + ndn = node.name + local_nodenames.append(ndn) + nlen = nlens[ndn] + + local_primal_residuals[ndn] = np.zeros(nlen, dtype="d") + global_primal_residuals[ndn] = np.zeros(nlen, dtype="d") + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + xbars = s._mpisppy_model.xbars + for node in s._mpisppy_node_list: + ndn = node.name + primal_residuals = local_primal_residuals[ndn] + + unweighted_primal_residuals = np.fromiter( + ( + abs(v._value - xbars[ndn, i]._value) + for i, v in enumerate(node.nonant_vardata_list) + ), + dtype="d", + count=nlens[ndn], + ) + primal_residuals += s._mpisppy_probability * unweighted_primal_residuals + + for nodename in local_nodenames: + ph.comms[nodename].Allreduce( + [local_primal_residuals[nodename], MPI.DOUBLE], + [global_primal_residuals[nodename], MPI.DOUBLE], + op=MPI.SUM, + ) + + primal_resid = {} + for ndn, global_primal_resid in global_primal_residuals.items(): + for i, v in enumerate(global_primal_resid): + primal_resid[ndn, i] = v + + return primal_resid + + @staticmethod + def _compute_min_max(ph, npop, mpiop, start): + local_nodenames = [] + local_xmaxmin = {} + global_xmaxmin = {} + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + for node in s._mpisppy_node_list: + if node.name not in local_nodenames: + ndn = node.name + local_nodenames.append(ndn) + nlen = nlens[ndn] + + local_xmaxmin[ndn] = start * np.ones(nlen, dtype="d") + global_xmaxmin[ndn] = np.zeros(nlen, dtype="d") + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + for node in s._mpisppy_node_list: + ndn = node.name + xmaxmin = local_xmaxmin[ndn] + + xmaxmin_partial = np.fromiter( + (v._value for v in node.nonant_vardata_list), + dtype="d", + count=nlens[ndn], + ) + xmaxmin = npop(xmaxmin, xmaxmin_partial) + local_xmaxmin[ndn] = xmaxmin + + for nodename in local_nodenames: + ph.comms[nodename].Allreduce( + [local_xmaxmin[nodename], MPI.DOUBLE], + [global_xmaxmin[nodename], MPI.DOUBLE], + op=mpiop, + ) + + xmaxmin_dict = {} + for ndn, global_xmaxmin_dict in global_xmaxmin.items(): + for i, v in enumerate(global_xmaxmin_dict): + xmaxmin_dict[ndn, i] = v + + return xmaxmin_dict + + @staticmethod + def _compute_xmax(ph): + return SepRho._compute_min_max(ph, np.maximum, MPI.MAX, -np.inf) + + @staticmethod + def _compute_xmin(ph): + return SepRho._compute_min_max(ph, np.minimum, MPI.MIN, np.inf) + + def pre_iter0(self): + pass + + def post_iter0(self): + ph = self.ph + primal_resid = self._compute_primal_residual_norm(ph) + xmax = self._compute_xmax(ph) + xmin = self._compute_xmin(ph) + + for s in ph.local_scenarios.values(): + cc = nonant_cost_coeffs(s) + for ndn_i, rho in s._mpisppy_model.rho.items(): + if cc[ndn_i] != 0: + nv = s._mpisppy_data.nonant_indices[ndn_i] # var_data object + if nv.is_integer(): + rho._value = abs(cc[ndn_i]) / (xmax[ndn_i] - xmin[ndn_i] + 1) + else: + rho._value = abs(cc[ndn_i]) / max(1, primal_resid[ndn_i]) + + rho._value *= self.multiplier + + # if ph.cylinder_rank==0: + # print(ndn_i,nv.getname(),xmax[ndn_i],xmin[ndn_i],primal_resid[ndn_i],cc[ndn_i],rho._value) + if ph.cylinder_rank == 0: + print("Rho values updated by SepRho Extension") + + def miditer(self): + pass + + def enditer(self): + pass + + def post_everything(self): + pass diff --git a/mpisppy/generic_cylinders.py b/mpisppy/generic_cylinders.py index 95200edbc..ebd009f24 100644 --- a/mpisppy/generic_cylinders.py +++ b/mpisppy/generic_cylinders.py @@ -70,6 +70,8 @@ def _parse_args(m): cfg.gradient_args() cfg.dynamic_gradient_args() cfg.reduced_costs_args() + cfg.sep_rho_args() + cfg.coeff_rho_args() cfg.parse_command_line(f"mpi-sppy for {cfg.module_name}") return cfg @@ -97,14 +99,17 @@ def _name_lists(module, cfg): def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement): rho_setter = module._rho_setter if hasattr(module, '_rho_setter') else None if cfg.default_rho is None and rho_setter is None: - raise RuntimeError("No rho_setter so a default must be specified via --default-rho") + if cfg.sep_rho or cfg.coeff_rho: + cfg.default_rho = 1 + else: + raise RuntimeError("No rho_setter so a default must be specified via --default-rho") - if False: # maybe later... cfg.use_norm_rho_converger: + if cfg.use_norm_rho_converger: if not cfg.use_norm_rho_updater: raise RuntimeError("--use-norm-rho-converger requires --use-norm-rho-updater") else: ph_converger = NormRhoConverger - elif False: # maybe later... cfg.primal_dual_converger: + elif cfg.primal_dual_converger: ph_converger = PrimalDualConverger else: ph_converger = None @@ -156,7 +161,13 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ if cfg.grad_rho_setter: ext_classes.append(Gradient_extension) hub_dict['opt_kwargs']['options']['gradient_extension_options'] = {'cfg': cfg} - + + if cfg.sep_rho: + vanilla.add_sep_rho(hub_dict, cfg) + + if cfg.coeff_rho: + vanilla.add_coeff_rho(hub_dict, cfg) + if len(ext_classes) != 0: hub_dict['opt_kwargs']['extensions'] = MultiExtension hub_dict["opt_kwargs"]["extension_kwargs"] = {"ext_classes" : ext_classes} @@ -179,6 +190,11 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ all_nodenames=all_nodenames, rho_setter=rho_setter, ) + # Need to fix FWPH to support extensions + # if cfg.sep_rho: + # vanilla.add_sep_rho(fw_spoke, cfg) + # if cfg.coeff_rho: + # vanilla.add_coeff_rho(fw_spoke, cfg) # Standard Lagrangian bound spoke if cfg.lagrangian: @@ -195,6 +211,10 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ rho_setter = rho_setter, all_nodenames = all_nodenames, ) + if cfg.sep_rho: + vanilla.add_sep_rho(ph_ob_spoke, cfg) + if cfg.coeff_rho: + vanilla.add_coeff_rho(ph_ob_spoke, cfg) # subgradient outer bound spoke if cfg.subgradient: @@ -203,6 +223,10 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ rho_setter = rho_setter, all_nodenames = all_nodenames, ) + if cfg.sep_rho: + vanilla.add_sep_rho(subgradient_spoke, cfg) + if cfg.coeff_rho: + vanilla.add_coeff_rho(subgradient_spoke, cfg) # xhat shuffle bound spoke if cfg.xhatshuffle: diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index 57ae43469..e09e8e675 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -40,6 +40,8 @@ from mpisppy.extensions.fixer import Fixer from mpisppy.extensions.cross_scen_extension import CrossScenarioExtension from mpisppy.extensions.reduced_costs_fixer import ReducedCostsFixer +from mpisppy.extensions.sep_rho import SepRho +from mpisppy.extensions.coeff_rho import CoeffRho from mpisppy.utils.wxbarreader import WXBarReader from mpisppy.utils.wxbarwriter import WXBarWriter @@ -105,6 +107,9 @@ def ph_hub( options["linearize_binary_proximal_terms"] = cfg.linearize_binary_proximal_terms options["linearize_proximal_terms"] = cfg.linearize_proximal_terms options["proximal_linearization_tolerance"] = cfg.proximal_linearization_tolerance + options["smoothed"] = 2 if cfg.smoothing else 0 + options["defaultPHp"] = cfg.smoothing_rho_ratio + options["defaultPHbeta"] = cfg.smoothing_beta hub_dict = { "hub_class": PHHub, @@ -196,6 +201,14 @@ def add_fixer(hub_dict, "id_fix_list_fct": cfg.id_fix_list_fct} return hub_dict +def add_sep_rho(hub_dict, cfg): + hub_dict = extension_adder(hub_dict,SepRho) + hub_dict["opt_kwargs"]["options"]["sep_rho_options"] = {"multiplier" : cfg.sep_rho_multiplier} + +def add_coeff_rho(hub_dict, cfg): + hub_dict = extension_adder(hub_dict,CoeffRho) + hub_dict["opt_kwargs"]["options"]["coeff_rho_options"] = {"multiplier" : cfg.coeff_rho_multiplier} + def add_cross_scenario_cuts(hub_dict, cfg, ): diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index 397b1179a..6955be408 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -248,6 +248,23 @@ def ph_args(self): domain=float, default=1.e-1) + self.add_to_config("smoothing", + description="For PH, add a smoothing term to the objective", + domain=bool, + default=False) + + self.add_to_config("smoothing_rho_ratio", + description="For PH, when smoothing, the ratio of " + "the smoothing coefficient to rho (default 1e-1)", + domain=float, + default=1.e-1) + + self.add_to_config("smoothing_beta", + description="For PH, when smoothing, the smoothing " + "memory coefficient beta (default 2e-1)", + domain=float, + default=2.e-1) + def make_parser(self, progname=None, num_scens_reqd=False): raise RuntimeError("make_parser is no longer used. See comments at top of config.py") @@ -393,6 +410,28 @@ def fixer_args(self): default=1e-2) + def sep_rho_args(self): + self.add_to_config("sep_rho", + description="have a SepRho extension", + domain=bool, + default=False) + self.add_to_config("sep_rho_multiplier", + description="multiplier for SepRho (default 1.0)", + domain=float, + default=1.0) + + + def coeff_rho_args(self): + self.add_to_config("coeff_rho", + description="have a CoeffRho extension", + domain=bool, + default=False) + self.add_to_config("coeff_rho_multiplier", + description="multiplier for CoeffRho (default 1.0)", + domain=float, + default=1.0) + + def gapper_args(self): self.add_to_config('mipgaps_json', diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py index b33639cf2..aa1866b47 100644 --- a/mpisppy/utils/sputils.py +++ b/mpisppy/utils/sputils.py @@ -16,6 +16,7 @@ import numpy as np import mpisppy.scenario_tree as scenario_tree from pyomo.core import Objective +from pyomo.repn import generate_standard_repn from mpisppy import MPI, haveMPI from pyomo.core.expr.numeric_expr import LinearExpression @@ -941,16 +942,52 @@ def reenable_tictoc_output(): tt_timer._ostream.close() tt_timer._ostream = sys.stdout - + def find_active_objective(pyomomodel): + return find_objective(pyomomodel, active=True) + + +def find_objective(pyomomodel, active=False): # return the only active objective or raise and error obj = list(pyomomodel.component_data_objects( Objective, active=True, descend_into=True)) - if len(obj) != 1: + if len(obj) == 1: + return obj[0] + if active or len(obj) > 1: raise RuntimeError("Could not identify exactly one active " "Objective for model '%s' (found %d objectives)" % (pyomomodel.name, len(obj))) - return obj[0] + # search again for a single inactive objective + obj = list(pyomomodel.component_data_objects( + Objective, descend_into=True)) + if len(obj) == 1: + return obj[0] + raise RuntimeError("Could not identify exactly one objective for model " + f"{pyomomodel.name} (found {len(obj)} objectives)") + + +def nonant_cost_coeffs(s): + """ + return a dictionary from s._mpisppy_data.nonant_indices.keys() + to the objective cost coefficient + """ + objective = find_objective(s) + + # initialize to 0 + cost_coefs = {ndn_i: 0 for ndn_i in s._mpisppy_data.nonant_indices} + repn = generate_standard_repn(objective.expr, quadratic=False) + for coef, var in zip(repn.linear_coefs, repn.linear_vars): + if id(var) in s._mpisppy_data.varid_to_nonant_index: + cost_coefs[s._mpisppy_data.varid_to_nonant_index[id(var)]] = coef + + for var in repn.nonlinear_vars: + if id(var) in s._mpisppy_data.varid_to_nonant_index: + raise RuntimeError( + "Found nonlinear variables in the objective function. " + f"Variable {var} has nonlinear interactions in the objective funtion" + ) + return cost_coefs + def create_nodenames_from_branching_factors(BFS): """