source: pyomo/trunk/examples/pysp/sizes/CCmodels/ReferenceModel.py @ 9239

Last change on this file since 9239 was 9239, checked in by wehart, 5 years ago

Converting files with lines containing coopr/Coopr -> pyomo/Pyomo

File size: 8.6 KB
Line 
1# revised to have a chance constraint for demand as needed by runef with the CC option
2#
3# This is the two-period version of the SIZES optimization model
4# derived from the three stage model in:
5#A. L{\o}kketangen and D. L. Woodruff,
6#"Progressive Hedging and Tabu Search Applied to Mixed Integer (0,1) Multistage Stochastic Programming",
7#Journal of Heuristics, 1996, Vol 2, Pages 111-128.
8
9from pyomo.core import *
10
11#
12# Model
13#
14
15model = AbstractModel()
16
17#
18# Parameters
19#
20
21# the number of product sizes.
22model.NumSizes = Param(within=NonNegativeIntegers)
23
24# the set of sizes, labeled 1 through NumSizes.
25def product_sizes_rule(model):
26    return set(range(1, model.NumSizes()+1))
27model.ProductSizes = Set(initialize=product_sizes_rule)
28
29# the deterministic demands for product at each size.
30model.DemandsFirstStage = Param(model.ProductSizes, within=NonNegativeIntegers)
31model.DemandsSecondStage = Param(model.ProductSizes, within=NonNegativeIntegers)
32
33# the unit production cost at each size.
34model.UnitProductionCosts = Param(model.ProductSizes, within=NonNegativeReals)
35
36# the setup cost for producing any units of size i.
37model.SetupCosts = Param(model.ProductSizes, within=NonNegativeReals)
38
39# the unit penalty cost of meeting demand for size j with larger size i.
40model.UnitPenaltyCosts = Param(model.ProductSizes, within=NonNegativeReals)
41
42# the cost to reduce a unit i to a lower unit j.
43model.UnitReductionCost = Param(within=NonNegativeReals)
44
45# a cap on the overall production within any time stage.
46model.Capacity = Param(within=PositiveReals)
47
48# a derived set to constrain the NumUnitsCut variable domain.
49# TBD: the (i,j) with i >= j set should be a generic utility.
50def num_units_cut_domain_rule(model):
51   ans = set()
52   for i in range(1,model.NumSizes()+1):
53      for j in range(1, i+1):
54         ans.add((i,j))   
55   return ans
56
57model.NumUnitsCutDomain = Set(initialize=num_units_cut_domain_rule, dimen=2)
58
59#
60# Variables
61#
62
63# are any products at size i produced?
64model.ProduceSizeFirstStage = Var(model.ProductSizes, domain=Boolean)
65model.ProduceSizeSecondStage = Var(model.ProductSizes, domain=Boolean)
66
67# NOTE: The following (num-produced and num-cut) variables are implicitly integer
68#       under the normal cost objective, but with the PH cost objective, this isn't
69#       the case.
70
71# the number of units at each size produced.
72model.NumProducedFirstStage = Var(model.ProductSizes, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
73model.NumProducedSecondStage = Var(model.ProductSizes, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
74
75# the number of units of size i cut (down) to meet demands for units of size j.
76model.NumUnitsCutFirstStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
77model.NumUnitsCutSecondStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity))
78
79# stage-specific cost variables for use in the pysp scenario tree / analysis.
80model.FirstStageCost = Var(domain=NonNegativeReals)
81model.SecondStageCost = Var(domain=NonNegativeReals)
82
83#
84# Constraints
85#
86
87# ensure that demand is satisfied in the first stage, accounting for cut-downs.
88def demand_satisfied_first_stage_rule(model, i):
89   return (0.0, sum([model.NumUnitsCutFirstStage[j,i] for j in model.ProductSizes if j >= i]) - model.DemandsFirstStage[i], None)
90model.DemandSatisfiedFirstStage = Constraint(model.ProductSizes, rule=demand_satisfied_first_stage_rule)
91
92###### CC: chance constraint #######
93## instead of requiring demand always be met in the second stage,
94## we set up a variable that indicates whether it is met
95##def demand_satisfied_second_stage_rule(i, model):
96##   return (0.0, sum([model.NumUnitsCutSecondStage[j,i] for j in model.ProductSizes if j >= i]) - model.DemandsSecondStage[i], None)   
97##model.DemandSatisfiedSecondStage = Constraint(model.ProductSizes, rule=demand_satisfied_second_stage_rule)
98
99model.dIndicator = Var(domain=Boolean)# indicate that all demans is met (for scenario)
100model.dforsize = Var(model.ProductSizes, domain=Boolean)  # indicate that demand is met for size
101
102# The production capacity per time stage serves as a simple upper bound for "M".
103def establish_dforsize_rule(model, i):
104   return model.dforsize[i]-1.0 <= (sum([model.NumUnitsCutSecondStage[j,i] for j in model.ProductSizes if j >= i]) - model.DemandsSecondStage[i]) / model.Capacity
105
106model.establish_dforsize = Constraint(model.ProductSizes, rule=establish_dforsize_rule)
107
108# it "wants to be a one" so don't let it unless it should be
109def establish_dIndicator_rule(model):
110   return model.dIndicator * model.NumSizes <= summation(model.dforsize)
111
112model.establish_dIndicator = Constraint(rule=establish_dIndicator_rule)
113
114# if the chance constraint is not desired, then runef can give the fully admissible solutions
115# by using the following AllDemandMet constraint
116# (this should be commented-out to get a chance constraint)
117#def AllDemandMet_rule(model):
118#   return model.dIndicator >= 1.0
119#model.AllDemandMet = Constraint(rule=AllDemandMet_rule)
120
121# ensure that you don't produce any units if the decision has been made to disable producion.
122def enforce_production_first_stage_rule(model, i):
123   # The production capacity per time stage serves as a simple upper bound for "M".
124   return (None, model.NumProducedFirstStage[i] - model.Capacity * model.ProduceSizeFirstStage[i], 0.0)
125
126def enforce_production_second_stage_rule(model, i):
127   # The production capacity per time stage serves as a simple upper bound for "M".   
128   return (None, model.NumProducedSecondStage[i] - model.Capacity * model.ProduceSizeSecondStage[i], 0.0)
129
130model.EnforceProductionBinaryFirstStage = Constraint(model.ProductSizes, rule=enforce_production_first_stage_rule)
131model.EnforceProductionBinarySecondStage = Constraint(model.ProductSizes, rule=enforce_production_second_stage_rule)
132
133# ensure that the production capacity is not exceeded for each time stage.
134def enforce_capacity_first_stage_rule(model):
135   return (None, sum([model.NumProducedFirstStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0)
136
137def enforce_capacity_second_stage_rule(model):
138   return (None, sum([model.NumProducedSecondStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0)   
139
140model.EnforceCapacityLimitFirstStage = Constraint(rule=enforce_capacity_first_stage_rule)
141model.EnforceCapacityLimitSecondStage = Constraint(rule=enforce_capacity_second_stage_rule)
142
143# ensure that you can't generate inventory out of thin air.
144def enforce_inventory_first_stage_rule(model, i):
145   return (None, \
146           sum([model.NumUnitsCutFirstStage[i,j] for j in model.ProductSizes if j <= i]) - \
147           model.NumProducedFirstStage[i], \
148           0.0)
149
150def enforce_inventory_second_stage_rule(model, i):
151   return (None, \
152           sum([model.NumUnitsCutFirstStage[i,j] for j in model.ProductSizes if j <= i]) + \
153           sum([model.NumUnitsCutSecondStage[i,j] for j in model.ProductSizes if j <= i]) \
154           - model.NumProducedFirstStage[i] - model.NumProducedSecondStage[i], \
155           0.0)
156
157model.EnforceInventoryFirstStage = Constraint(model.ProductSizes, rule=enforce_inventory_first_stage_rule)
158model.EnforceInventorySecondStage = Constraint(model.ProductSizes, rule=enforce_inventory_second_stage_rule)
159
160# stage-specific cost computations.
161def first_stage_cost_rule(model):
162   production_costs = sum([model.SetupCosts[i] * model.ProduceSizeFirstStage[i] + \
163                          model.UnitProductionCosts[i] * model.NumProducedFirstStage[i] \
164                          for i in model.ProductSizes])
165   cut_costs = sum([model.UnitReductionCost * model.NumUnitsCutFirstStage[i,j] \
166                   for (i,j) in model.NumUnitsCutDomain if i != j])
167   return (model.FirstStageCost - production_costs - cut_costs) == 0.0
168
169model.ComputeFirstStageCost = Constraint(rule=first_stage_cost_rule)
170
171def second_stage_cost_rule(model):
172   production_costs = sum([model.SetupCosts[i] * model.ProduceSizeSecondStage[i] + \
173                          model.UnitProductionCosts[i] * model.NumProducedSecondStage[i] \
174                          for i in model.ProductSizes])
175   cut_costs = sum([model.UnitReductionCost * model.NumUnitsCutSecondStage[i,j] \
176                   for (i,j) in model.NumUnitsCutDomain if i != j])
177   return (model.SecondStageCost - production_costs - cut_costs) == 0.0   
178
179model.ComputeSecondStageCost = Constraint(rule=second_stage_cost_rule)
180
181#
182# PySP Auto-generated Objective
183#
184# minimize: sum of StageCostVariables
185#
186# A active scenario objective equivalent to that generated by PySP is
187# included here for informational purposes.
188def total_cost_rule(model):
189    return model.FirstStageCost + model.SecondStageCost
190model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)
191
Note: See TracBrowser for help on using the repository browser.