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

Last change on this file since 9501 was 9501, checked in by jwatson, 4 years ago

Updated chance-constrained sizes model.

File size: 8.7 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
101model.lambdaMult = Param(initialize=0.0, mutable=True)
102
103# The production capacity per time stage serves as a simple upper bound for "M".
104def establish_dforsize_rule(model, i):
105   return model.dforsize[i]-1.0 <= (sum([model.NumUnitsCutSecondStage[j,i] for j in model.ProductSizes if j >= i]) - model.DemandsSecondStage[i]) / model.Capacity
106
107model.establish_dforsize = Constraint(model.ProductSizes, rule=establish_dforsize_rule)
108
109# it "wants to be a one" so don't let it unless it should be
110def establish_dIndicator_rule(model):
111   return model.dIndicator * model.NumSizes <= summation(model.dforsize)
112
113model.establish_dIndicator = Constraint(rule=establish_dIndicator_rule)
114
115# if the chance constraint is not desired, then runef can give the fully admissible solutions
116# by using the following AllDemandMet constraint
117# (this should be commented-out to get a chance constraint)
118##def AllDemandMet_rule(model):
119##   return model.dIndicator >= 1.0
120##model.AllDemandMet = Constraint(rule=AllDemandMet_rule)
121
122# ensure that you don't produce any units if the decision has been made to disable producion.
123def enforce_production_first_stage_rule(model, i):
124   # The production capacity per time stage serves as a simple upper bound for "M".
125   return (None, model.NumProducedFirstStage[i] - model.Capacity * model.ProduceSizeFirstStage[i], 0.0)
126
127def enforce_production_second_stage_rule(model, i):
128   # The production capacity per time stage serves as a simple upper bound for "M".   
129   return (None, model.NumProducedSecondStage[i] - model.Capacity * model.ProduceSizeSecondStage[i], 0.0)
130
131model.EnforceProductionBinaryFirstStage = Constraint(model.ProductSizes, rule=enforce_production_first_stage_rule)
132model.EnforceProductionBinarySecondStage = Constraint(model.ProductSizes, rule=enforce_production_second_stage_rule)
133
134# ensure that the production capacity is not exceeded for each time stage.
135def enforce_capacity_first_stage_rule(model):
136   return (None, sum([model.NumProducedFirstStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0)
137
138def enforce_capacity_second_stage_rule(model):
139   return (None, sum([model.NumProducedSecondStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0)   
140
141model.EnforceCapacityLimitFirstStage = Constraint(rule=enforce_capacity_first_stage_rule)
142model.EnforceCapacityLimitSecondStage = Constraint(rule=enforce_capacity_second_stage_rule)
143
144# ensure that you can't generate inventory out of thin air.
145def enforce_inventory_first_stage_rule(model, i):
146   return (None, \
147           sum([model.NumUnitsCutFirstStage[i,j] for j in model.ProductSizes if j <= i]) - \
148           model.NumProducedFirstStage[i], \
149           0.0)
150
151def enforce_inventory_second_stage_rule(model, i):
152   return (None, \
153           sum([model.NumUnitsCutFirstStage[i,j] for j in model.ProductSizes if j <= i]) + \
154           sum([model.NumUnitsCutSecondStage[i,j] for j in model.ProductSizes if j <= i]) \
155           - model.NumProducedFirstStage[i] - model.NumProducedSecondStage[i], \
156           0.0)
157
158model.EnforceInventoryFirstStage = Constraint(model.ProductSizes, rule=enforce_inventory_first_stage_rule)
159model.EnforceInventorySecondStage = Constraint(model.ProductSizes, rule=enforce_inventory_second_stage_rule)
160
161# stage-specific cost computations.
162def first_stage_cost_rule(model):
163   production_costs = sum([model.SetupCosts[i] * model.ProduceSizeFirstStage[i] + \
164                          model.UnitProductionCosts[i] * model.NumProducedFirstStage[i] \
165                          for i in model.ProductSizes])
166   cut_costs = sum([model.UnitReductionCost * model.NumUnitsCutFirstStage[i,j] \
167                   for (i,j) in model.NumUnitsCutDomain if i != j])
168   return (model.FirstStageCost - production_costs - cut_costs) == 0.0
169
170model.ComputeFirstStageCost = Constraint(rule=first_stage_cost_rule)
171
172def second_stage_cost_rule(model):
173   production_costs = sum([model.SetupCosts[i] * model.ProduceSizeSecondStage[i] + \
174                          model.UnitProductionCosts[i] * model.NumProducedSecondStage[i] \
175                          for i in model.ProductSizes])
176   cut_costs = sum([model.UnitReductionCost * model.NumUnitsCutSecondStage[i,j] \
177                   for (i,j) in model.NumUnitsCutDomain if i != j])
178   return (model.SecondStageCost - production_costs - cut_costs) == 0.0   
179
180model.ComputeSecondStageCost = Constraint(rule=second_stage_cost_rule)
181
182#
183# PySP Auto-generated Objective
184#
185# minimize: sum of StageCostVariables
186#
187# A active scenario objective equivalent to that generated by PySP is
188# included here for informational purposes.
189def total_cost_rule(model):
190    return model.FirstStageCost + model.SecondStageCost
191model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)
192
Note: See TracBrowser for help on using the repository browser.