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 | |
---|
9 | from pyomo.core import * |
---|
10 | |
---|
11 | # |
---|
12 | # Model |
---|
13 | # |
---|
14 | |
---|
15 | model = AbstractModel() |
---|
16 | |
---|
17 | # |
---|
18 | # Parameters |
---|
19 | # |
---|
20 | |
---|
21 | # the number of product sizes. |
---|
22 | model.NumSizes = Param(within=NonNegativeIntegers) |
---|
23 | |
---|
24 | # the set of sizes, labeled 1 through NumSizes. |
---|
25 | def product_sizes_rule(model): |
---|
26 | return set(range(1, model.NumSizes()+1)) |
---|
27 | model.ProductSizes = Set(initialize=product_sizes_rule) |
---|
28 | |
---|
29 | # the deterministic demands for product at each size. |
---|
30 | model.DemandsFirstStage = Param(model.ProductSizes, within=NonNegativeIntegers) |
---|
31 | model.DemandsSecondStage = Param(model.ProductSizes, within=NonNegativeIntegers) |
---|
32 | |
---|
33 | # the unit production cost at each size. |
---|
34 | model.UnitProductionCosts = Param(model.ProductSizes, within=NonNegativeReals) |
---|
35 | |
---|
36 | # the setup cost for producing any units of size i. |
---|
37 | model.SetupCosts = Param(model.ProductSizes, within=NonNegativeReals) |
---|
38 | |
---|
39 | # the unit penalty cost of meeting demand for size j with larger size i. |
---|
40 | model.UnitPenaltyCosts = Param(model.ProductSizes, within=NonNegativeReals) |
---|
41 | |
---|
42 | # the cost to reduce a unit i to a lower unit j. |
---|
43 | model.UnitReductionCost = Param(within=NonNegativeReals) |
---|
44 | |
---|
45 | # a cap on the overall production within any time stage. |
---|
46 | model.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. |
---|
50 | def 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 | |
---|
57 | model.NumUnitsCutDomain = Set(initialize=num_units_cut_domain_rule, dimen=2) |
---|
58 | |
---|
59 | # |
---|
60 | # Variables |
---|
61 | # |
---|
62 | |
---|
63 | # are any products at size i produced? |
---|
64 | model.ProduceSizeFirstStage = Var(model.ProductSizes, domain=Boolean) |
---|
65 | model.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. |
---|
72 | model.NumProducedFirstStage = Var(model.ProductSizes, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity)) |
---|
73 | model.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. |
---|
76 | model.NumUnitsCutFirstStage = Var(model.NumUnitsCutDomain, domain=NonNegativeIntegers, bounds=(0.0, model.Capacity)) |
---|
77 | model.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. |
---|
80 | model.FirstStageCost = Var(domain=NonNegativeReals) |
---|
81 | model.SecondStageCost = Var(domain=NonNegativeReals) |
---|
82 | |
---|
83 | # |
---|
84 | # Constraints |
---|
85 | # |
---|
86 | |
---|
87 | # ensure that demand is satisfied in the first stage, accounting for cut-downs. |
---|
88 | def 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) |
---|
90 | model.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 | |
---|
99 | model.dIndicator = Var(domain=Boolean)# indicate that all demans is met (for scenario) |
---|
100 | model.dforsize = Var(model.ProductSizes, domain=Boolean) # indicate that demand is met for size |
---|
101 | model.lambdaMult = Param(initialize=0.0, mutable=True) |
---|
102 | |
---|
103 | # The production capacity per time stage serves as a simple upper bound for "M". |
---|
104 | def 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 | |
---|
107 | model.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 |
---|
110 | def establish_dIndicator_rule(model): |
---|
111 | return model.dIndicator * model.NumSizes <= summation(model.dforsize) |
---|
112 | |
---|
113 | model.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. |
---|
123 | def 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 | |
---|
127 | def 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 | |
---|
131 | model.EnforceProductionBinaryFirstStage = Constraint(model.ProductSizes, rule=enforce_production_first_stage_rule) |
---|
132 | model.EnforceProductionBinarySecondStage = Constraint(model.ProductSizes, rule=enforce_production_second_stage_rule) |
---|
133 | |
---|
134 | # ensure that the production capacity is not exceeded for each time stage. |
---|
135 | def enforce_capacity_first_stage_rule(model): |
---|
136 | return (None, sum([model.NumProducedFirstStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0) |
---|
137 | |
---|
138 | def enforce_capacity_second_stage_rule(model): |
---|
139 | return (None, sum([model.NumProducedSecondStage[i] for i in model.ProductSizes]) - model.Capacity, 0.0) |
---|
140 | |
---|
141 | model.EnforceCapacityLimitFirstStage = Constraint(rule=enforce_capacity_first_stage_rule) |
---|
142 | model.EnforceCapacityLimitSecondStage = Constraint(rule=enforce_capacity_second_stage_rule) |
---|
143 | |
---|
144 | # ensure that you can't generate inventory out of thin air. |
---|
145 | def 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 | |
---|
151 | def 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 | |
---|
158 | model.EnforceInventoryFirstStage = Constraint(model.ProductSizes, rule=enforce_inventory_first_stage_rule) |
---|
159 | model.EnforceInventorySecondStage = Constraint(model.ProductSizes, rule=enforce_inventory_second_stage_rule) |
---|
160 | |
---|
161 | # stage-specific cost computations. |
---|
162 | def 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 | |
---|
170 | model.ComputeFirstStageCost = Constraint(rule=first_stage_cost_rule) |
---|
171 | |
---|
172 | def 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 | |
---|
180 | model.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. |
---|
189 | def total_cost_rule(model): |
---|
190 | return model.FirstStageCost + model.SecondStageCost |
---|
191 | model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize) |
---|
192 | |
---|