Merge pull request #483 from DLWoodruff/sizesEval

.ruff.toml

@@ -4,4 +4,5 @@ extend-exclude = [
     "./mpisppy/utils/callbacks/termination/tests/markshare2.py",
     "./mpisppy/tests/examples/hydro/hydro.py",
     "./examples/hydro/hydro.py",
+    "./examples/sizes/models/ExpressionModel.py",
 ]

examples/sizes/models/ExpressionModel.py

@@ -0,0 +1,182 @@
+###############################################################################
+# 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.
+###############################################################################
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and 
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain 
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+#
+# This is the two-period version of the SIZES optimization model
+# derived from the three stage model in:
+#A. L{\o}kketangen and D. L. Woodruff,
+#"Progressive Hedging and Tabu Search Applied to Mixed Integer (0,1) Multistage Stochastic Programming",
+#Journal of Heuristics, 1996, Vol 2, Pages 111-128.
+# This version uses expressions for stage costs. January 2025.
+
+from pyomo.core import *
+
+#
+# Model
+#
+
+model = AbstractModel()
+
+#
+# Parameters
+#
+
+# the number of product sizes.
+model.NumSizes = Param(within=NonNegativeIntegers)
+
+# the set of sizes, labeled 1 through NumSizes.
+def product_sizes_rule(model):
+    return list(range(1, model.NumSizes()+1))
+model.ProductSizes = Set(initialize=product_sizes_rule)
+
+# the deterministic demands for product at each size.
+model.DemandsFirstStage = Param(model.ProductSizes, within=NonNegativeIntegers)
+model.DemandsSecondStage = Param(model.ProductSizes, within=NonNegativeIntegers)
+
+# the unit production cost at each size.
+model.UnitProductionCosts = Param(model.ProductSizes, within=NonNegativeReals)
+
+# the setup cost for producing any units of size i.
+model.SetupCosts = Param(model.ProductSizes, within=NonNegativeReals)
+
+# the cost to reduce a unit i to a lower unit j.
+model.UnitReductionCost = Param(within=NonNegativeReals)
+
+# a cap on the overall production within any time stage.
+model.Capacity = Param(within=PositiveReals)
+
+# a derived set to constrain the NumUnitsCut variable domain.
+# TBD: the (i,j) with i >= j set should be a generic utility.
+def num_units_cut_domain_rule(model):
+    ans = list()
+    for i in range(1,model.NumSizes()+1):
+        for j in range(1, i+1):
+            ans.append((i,j))
+    return ans
+
+model.NumUnitsCutDomain = Set(initialize=num_units_cut_domain_rule, dimen=2)
+
+#
+# Variables
+#
+
+# are any products at size i produced?
+model.ProduceSizeFirstStage = Var(model.ProductSizes, domain=Boolean)
+model.ProduceSizeSecondStage = Var(model.ProductSizes, domain=Boolean)
+
+# NOTE: The following (num-produced and num-cut) variables are implicitly integer
+#       under the normal cost objective, but with the PH cost objective, this isn't
+#       the case.
+
+# the number of units at each size produced.
+model.NumProducedFirstStage = Var(model.ProductSizes, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
+model.NumProducedSecondStage = Var(model.ProductSizes, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
+
+# the number of units of size i cut (down) to meet demand for units of size j.
+model.NumUnitsCutFirstStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
+model.NumUnitsCutSecondStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
+
+# stage-specific cost variables for use in the pysp scenario tree / analysis.
+# changed to Expressions January 2025
+#model.FirstStageCost = Var(domain=NonNegativeReals)
+#model.SecondStageCost = Var(domain=NonNegativeReals)
+
+#
+# Constraints
+#
+
+# ensure that demand is satisfied in each time stage, accounting for cut-downs.
+def demand_satisfied_first_stage_rule(model, i):
+    return (0.0, sum([model.NumUnitsCutFirstStage[j,i] for j in model.ProductSizes if j >= i]) - model.DemandsFirstStage[i], None)
+
+def demand_satisfied_second_stage_rule(model, i):
+    return (0.0, sum([model.NumUnitsCutSecondStage[j,i] for j in model.ProductSizes if j >= i]) - model.DemandsSecondStage[i], None)
+
+model.DemandSatisfiedFirstStage = Constraint(model.ProductSizes, rule=demand_satisfied_first_stage_rule)
+model.DemandSatisfiedSecondStage = Constraint(model.ProductSizes, rule=demand_satisfied_second_stage_rule)
+
+# ensure that you don't produce any units if the decision has been made to disable producion.
+def enforce_production_first_stage_rule(model, i):
+    # The production capacity per time stage serves as a simple upper bound for "M".
+    return (None, model.NumProducedFirstStage[i] - model.Capacity * model.ProduceSizeFirstStage[i], 0.0)
+
+def enforce_production_second_stage_rule(model, i):
+    # The production capacity per time stage serves as a simple upper bound for "M".
+    return (None, model.NumProducedSecondStage[i] - model.Capacity * model.ProduceSizeSecondStage[i], 0.0)
+
+model.EnforceProductionBinaryFirstStage = Constraint(model.ProductSizes, rule=enforce_production_first_stage_rule)
+model.EnforceProductionBinarySecondStage = Constraint(model.ProductSizes, rule=enforce_production_second_stage_rule)
+
+# ensure that the production capacity is not exceeded for each time stage.
+def enforce_capacity_first_stage_rule(model):
+    return (None, sum([model.NumProducedFirstStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0)
+
+def enforce_capacity_second_stage_rule(model):
+    return (None, sum([model.NumProducedSecondStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0)
+
+model.EnforceCapacityLimitFirstStage = Constraint(rule=enforce_capacity_first_stage_rule)
+model.EnforceCapacityLimitSecondStage = Constraint(rule=enforce_capacity_second_stage_rule)
+
+# ensure that you can't generate inventory out of thin air.
+def enforce_inventory_first_stage_rule(model, i):
+    return (None, \
+            sum([model.NumUnitsCutFirstStage[i,j] for j in model.ProductSizes if j <= i]) - \
+            model.NumProducedFirstStage[i], \
+            0.0)
+
+def enforce_inventory_second_stage_rule(model, i):
+    return (None, \
+            sum([model.NumUnitsCutFirstStage[i,j] for j in model.ProductSizes if j <= i]) + \
+            sum([model.NumUnitsCutSecondStage[i,j] for j in model.ProductSizes if j <= i]) \
+            - model.NumProducedFirstStage[i] - model.NumProducedSecondStage[i], \
+            0.0)
+
+model.EnforceInventoryFirstStage = Constraint(model.ProductSizes, rule=enforce_inventory_first_stage_rule)
+model.EnforceInventorySecondStage = Constraint(model.ProductSizes, rule=enforce_inventory_second_stage_rule)
+
+# stage-specific cost computations.
+def first_stage_cost_rule(model):
+    production_costs = sum([model.SetupCosts[i] * model.ProduceSizeFirstStage[i] + \
+                           model.UnitProductionCosts[i] * model.NumProducedFirstStage[i] \
+                           for i in model.ProductSizes])
+    cut_costs = sum([model.UnitReductionCost * model.NumUnitsCutFirstStage[i,j] \
+                    for (i,j) in model.NumUnitsCutDomain if i != j])
+    return (production_costs - cut_costs)
+
+model.FirstStageCost = Expression(rule=first_stage_cost_rule)
+
+def second_stage_cost_rule(model):
+    production_costs = sum([model.SetupCosts[i] * model.ProduceSizeSecondStage[i] + \
+                           model.UnitProductionCosts[i] * model.NumProducedSecondStage[i] \
+                           for i in model.ProductSizes])
+    cut_costs = sum([model.UnitReductionCost * model.NumUnitsCutSecondStage[i,j] \
+                    for (i,j) in model.NumUnitsCutDomain if i != j])
+    return (production_costs - cut_costs)
+
+model.SecondStageCost = Expression(rule=second_stage_cost_rule)
+
+#
+# PySP Auto-generated Objective
+#
+# minimize: sum of StageCosts
+#
+# A active scenario objective equivalent to that generated by PySP is
+# included here for informational purposes.
+def total_cost_rule(model):
+    return model.FirstStageCost + model.SecondStageCost
+model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)
+

examples/sizes/sizes_expression.py

@@ -0,0 +1,164 @@
+###############################################################################
+# 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 os
+import models.ExpressionModel as ref  # the version with stage cost expressions
+import mpisppy.utils.sputils as sputils
+
+def scenario_creator(scenario_name, scenario_count=None):
+    if scenario_count not in (3, 10):
+        raise RuntimeError(
+            "scenario_count passed to scenario counter must equal either 3 or 10"
+        )
+
+    sizes_dir = os.path.dirname(__file__)
+    datadir = os.sep.join((sizes_dir, f"SIZES{scenario_count}"))
+    try:
+        fname = datadir + os.sep + scenario_name + ".dat"
+    except Exception:
+        print("FAIL: datadir=", datadir, " scenario_name=", scenario_name)
+
+    model = ref.model.create_instance(fname)
+
+    # now attach the one and only tree node
+    varlist = [model.NumProducedFirstStage, model.NumUnitsCutFirstStage]
+    sputils.attach_root_node(model, model.FirstStageCost, varlist)
+
+    return model
+
+
+def scenario_denouement(rank, scenario_name, scenario):
+    pass
+
+########## helper functions ########
+
+#=========
+def scenario_names_creator(num_scens,start=None):
+    # if start!=None, the list starts with the 'start' labeled scenario
+    # note that the scenarios for the sizes problem are one-based
+    if (start is None) :
+        start=1
+    return [f"Scenario{i}" for i in range(start, start+num_scens)]
+
+
+#=========
+def inparser_adder(cfg):
+    # add options unique to sizes
+    cfg.num_scens_required()
+    cfg.mip_options()
+
+
+#=========
+def kw_creator(cfg):
+    # (for Amalgamator): linked to the scenario_creator and inparser_adder
+    if cfg.num_scens not in (3, 10):
+        raise RuntimeError(f"num_scen must the 3 or 10; was {cfg.num_scen}")    
+    kwargs = {"scenario_count": cfg.num_scens}
+    return kwargs
+
+def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed,
+                             given_scenario=None, **scenario_creator_kwargs):
+    """ Create a scenario within a sample tree. Mainly for multi-stage and simple for two-stage.
+        (this function supports zhat and confidence interval code)
+    Args:
+        sname (string): scenario name to be created
+        stage (int >=1 ): for stages > 1, fix data based on sname in earlier stages
+        sample_branching_factors (list of ints): branching factors for the sample tree
+        seed (int): To allow random sampling (for some problems, it might be scenario offset)
+        given_scenario (Pyomo concrete model): if not None, use this to get data for ealier stages
+        scenario_creator_kwargs (dict): keyword args for the standard scenario creator funcion
+    Returns:
+        scenario (Pyomo concrete model): A scenario for sname with data in stages < stage determined
+                                         by the arguments
+    """
+    # Since this is a two-stage problem, we don't have to do much.
+    sca = scenario_creator_kwargs.copy()
+    sca["seedoffset"] = seed
+    sca["num_scens"] = sample_branching_factors[0]  # two-stage problem
+    return scenario_creator(sname, **sca)
+
+######## end helper functions #########
+
+########## a customized rho setter #############
+# If you are using sizes.py as a starting point for your model,
+#  you should be aware that you don't need a _rho_setter function.
+# This demonstrates how to use a customized rho setter; consider instead
+#  a gradient based rho setter.
+# note that _rho_setter is a reserved name....
+
+def _rho_setter(scen, **kwargs):
+    """ rho values for the scenario.
+    Args:
+        scen (pyo.ConcreteModel): the scenario
+    Returns:
+        a list of (id(vardata), rho)
+    Note:
+        This rho_setter will not work with proper bundles.
+    """
+    retlist = []
+    if not hasattr(scen, "UnitReductionCost"):
+        print("WARNING: _rho_setter not used (probably because of proper bundles)")
+        return retlist        
+    RF = 0.001  # a factor for rho, if you like
+
+    if "RF" in kwargs and isinstance(kwargs["RF"], float):
+        RF = kwargs["RF"]
+
+    cutrho = scen.UnitReductionCost * RF
+
+    for i in scen.ProductSizes:
+        idv = id(scen.NumProducedFirstStage[i])
+        rho = scen.UnitProductionCosts[i] * RF
+        retlist.append((idv, rho))
+
+        for j in scen.ProductSizes:
+            if j <= i:
+                idv = id(scen.NumUnitsCutFirstStage[i, j])
+                retlist.append((idv, cutrho))
+
+    return retlist
+
+
+def id_fix_list_fct(s):
+    """ specify tuples used by the fixer.
+
+    Args:
+        s (ConcreteModel): the sizes instance.
+    Returns:
+         i0, ik (tuples): one for iter 0 and other for general iterations.
+             Var id,  threshold, nb, lb, ub
+             The threshold is on the square root of the xbar squared differnce
+             nb, lb an bu an "no bound", "upper" and "lower" and give the numver
+                 of iterations or None for ik and for i0 anything other than None
+                 or None. In both cases, None indicates don't fix.
+    """
+    import mpisppy.extensions.fixer as fixer
+
+    iter0tuples = []
+    iterktuples = []
+    for i in s.ProductSizes:
+        iter0tuples.append(
+            fixer.Fixer_tuple(s.NumProducedFirstStage[i], th=0.01, nb=None, lb=0, ub=0)
+        )
+        iterktuples.append(
+            fixer.Fixer_tuple(s.NumProducedFirstStage[i], th=0.2, nb=3, lb=1, ub=2)
+        )
+        for j in s.ProductSizes:
+            if j <= i:
+                iter0tuples.append(
+                    fixer.Fixer_tuple(
+                        s.NumUnitsCutFirstStage[i, j], th=0.5, nb=None, lb=0, ub=0
+                    )
+                )
+                iterktuples.append(
+                    fixer.Fixer_tuple(
+                        s.NumUnitsCutFirstStage[i, j], th=0.2, nb=3, lb=1, ub=2
+                    )
+                )
+
+    return iter0tuples, iterktuples