source: trunk/examples/pysp/forestry/models-nb-yr/ReferenceModel.py @ 1759

Last change on this file since 1759 was 1759, checked in by jwatson, 10 years ago

Re-structuring the naming convention for the PySP forestry example.

  • Property svn:executable set to *
File size: 7.7 KB
Line 
1# fbv - Chile rules!
2# dlw - base forest model - begin April 2009
3# mnr - base forest model - revised May 1, 2009 - changed the notation of parameters to match those in paper.
4
5#
6# Imports (i.e., magical incantations)
7#
8
9import sys
10import os
11from os.path import abspath, dirname
12sys.path.insert(0, dirname(dirname(dirname(dirname(abspath(__file__))))))
13from coopr.pyomo import *
14
15#
16# Model
17#
18model = Model()
19model.name = "Forest Management Base Model"
20
21#
22# Implementation Sets (not in paper, used to make it happen)
23# general nodes (not in paper)
24model.Nodes = Set()
25
26#
27# Sets
28#
29
30# time Horizon \mathcal{T}
31model.Times = Set(ordered=True)
32
33# network nodes \mathcal{O}
34model.OriginNodes = Set(within=model.Nodes)
35
36# intersection nodes
37model.IntersectionNodes = Set(within=model.Nodes)
38
39# wood exit nodes
40model.ExitNodes = Set(within=model.Nodes)
41
42# harvest cells \mathcal{H}
43model.HarvestCells = Set()
44
45# harvest cells for origin o \mathcal{H}_o
46model.HCellsForOrigin = Set(model.OriginNodes, within=model.HarvestCells)
47
48# existing roads
49model.ExistingRoads = Set(within=model.Nodes*model.Nodes) 
50
51# potential roads
52model.PotentialRoads = Set(within=model.Nodes*model.Nodes)
53
54# Roads declaring
55#model.AllRoads = model.ExistingRoads & model.PotentialRoads
56model.AllRoads = Set(within=model.Nodes*model.Nodes)
57
58# derived set
59def aplus_arc_set_rule(v, model):
60   return  [(i, j) for (i, j) in model.AllRoads if j == v]
61model.Aplus = Set(model.Nodes, within=model.AllRoads, rule=aplus_arc_set_rule)
62
63# derived set
64def aminus_arc_set_rule(v, model):
65   return  [(i, j) for (i, j) in model.AllRoads if i == v]
66model.Aminus = Set(model.Nodes, within=model.AllRoads, rule=aminus_arc_set_rule)
67
68#
69# Deterministic Parameters
70#
71
72# productivity of cell h if it is harvested in period t
73model.a = Param(model.HarvestCells, model.Times)
74
75# Area of cell h to be harvested
76model.A = Param(model.HarvestCells)
77
78# flow capacity (smallest big M)
79def Umax_init(model):
80   return (sum([1.1 * model.a[h, t] * model.A[h] for h in model.HarvestCells for t in model.Times])/4)()
81model.Umax = Param(initialize=Umax_init)
82
83# harvesting cost of one hectare of cell h in time t
84model.P = Param(model.HarvestCells, model.Times)
85
86# unit Production cost at origin o in time t
87model.Q = Param(model.OriginNodes, model.Times)
88
89# construction cost of one road in arc (k,l) in time t
90model.C = Param(model.PotentialRoads, model.Times, default=0.0)
91
92# unit transport cost through arc (k,l) in time t
93model.D = Param(model.AllRoads, model.Times)
94
95#
96# Stochastic Parameters
97#
98
99# sale price at exit s in time period t
100model.R = Param(model.ExitNodes, model.Times)
101
102# upper and lower bounds on wood supplied
103model.Zlb = Param(model.Times)
104model.Zub = Param(model.Times)
105
106# yield ratio of the forest
107model.yr = Param(model.Times)
108
109#
110# Variables
111#
112
113# If cell h is harvested in period t
114model.delta = Var(model.HarvestCells, model.Times, domain=Boolean)
115
116# If road in arc (k,l) is built in period t
117model.gamma = Var(model.PotentialRoads, model.Times, domain=Boolean)
118
119# Flow of wood thransported through arc (k,l)
120model.f = Var(model.AllRoads, model.Times, domain=NonNegativeReals)
121
122# Supply of wood at exit s in period t
123model.z = Var(model.ExitNodes, model.Times, domain=NonNegativeReals)
124
125# Declared when changing the program for stochastics
126model.AnoProfit = Var(model.Times)
127
128#
129# Constraints
130#
131
132# flow balance at origin nodes
133def origin_flow_bal(o, t, model):
134   ans1 = sum([model.a[h, t] * model.yr[t] * model.A[h] * model.delta[h, t] for h in model.HCellsForOrigin[o]])
135   ans2 = sum([model.f[k, o, t] for (k, o) in model.Aplus[o]])
136   ans3 = sum([model.f[o, k, t] for (o, k) in model.Aminus[o]])
137   return (ans1 + ans2 - ans3) == 0.0
138
139model.EnforceOriginFlowBalance = Constraint(model.OriginNodes, model.Times, rule=origin_flow_bal)
140
141# flow balance at intersection nodes
142def intersection_flow_bal(j, t, model):
143   return (sum([model.f[k, j, t] for (k, j) in model.Aplus[j]]) - sum([model.f[j, k, t] for (j, k) in model.Aminus[j]])) == 0.0
144
145model.EnforceIntersectionFlowBalance = Constraint(model.IntersectionNodes, model.Times, rule=intersection_flow_bal)
146
147# flow balance at destination nodes
148def destination_flow_bal(e, t, model):
149   return (model.z[e, t] - sum([model.f[k, e, t] for (k, e) in model.Aplus[e]]) + sum([model.f[e, k, t] for (e, k) in model.Aminus[e]])) == 0.0
150
151model.EnforceDestinationFlowBalance = Constraint(model.ExitNodes, model.Times, rule=destination_flow_bal)
152
153# explicit divergence or supply equals demand or no inventory between periods
154def divergence_flow_bal(t, model):
155   return (sum([model.z[e, t] for e in model.ExitNodes]) - sum([model.a[h, t] * model.yr[t] * model.A[h] * model.delta[h, t] for h in model.HarvestCells])) == 0.0
156
157model.EnforceDivergenceFlowBalance = Constraint(model.Times, rule=divergence_flow_bal)
158
159# Wood Production Bounds
160def wood_prod_bounds(t, model):
161   return (model.Zlb[t], sum([model.z[e, t] for e in model.ExitNodes]), model.Zub[t])
162
163model.EnforceWoodProductionBounds = Constraint(model.Times, rule=wood_prod_bounds)
164
165# potential roads flow capacity 1
166def pot_road_flow_capacity(k, l, t, model):
167        return (None, model.f[k, l, t] - sum([model.Umax * model.gamma[k, l, tau] for tau in model.Times if tau <= t]), 0.0)
168
169model.EnforcePotentialRoadFlowCap = Constraint(model.PotentialRoads, model.Times, rule=pot_road_flow_capacity)
170
171# each potential road can be built no more than once in the time horizon
172def one_potential_road(k, l, model):
173   return (None, sum([model.gamma[k ,l, t] for t in model.Times]), 1.0)
174
175model.EnforceOnePotentialRoad = Constraint(model.PotentialRoads, rule=one_potential_road)
176
177# existing roads flow capacity
178def existing_roads_capacity(k, l, t, model):
179        return (None, model.f[k, l, t], model.Umax)
180
181model.EnforceExistingRoadCapacities = Constraint(model.ExistingRoads, model.Times, rule=existing_roads_capacity)
182
183# additional road bound capacity
184def additional_roads_capacity(k, l, t, model):
185        return ( None, model.f[k, l, t] - sum([model.z[e, t] for e in model.ExitNodes]), 0.0 )
186
187model.EnforceAdditionalRoadCapacities = Constraint(model.AllRoads, model.Times, rule=additional_roads_capacity)
188
189# each cell can be harvested no more than once in the time horizon
190def one_harvest(h, model):
191   return (None, sum([model.delta[h,t] for t in model.Times]), 1.0)
192
193model.EnforceOneHarvest = Constraint(model.HarvestCells, rule=one_harvest)
194
195#
196# Stage-specific profit computations
197#
198
199def stage_profit_rule(t, model):
200    # ans1 is T1 (i.e. wood sale revenue)
201    ans1 = sum([model.R[e, t] * model.z[e, t] for e in model.ExitNodes])
202   
203    # ans2 is T2 (i.e. wood harvest cost)
204    ans2 = sum([model.P[h, t] * model.A[h] * model.delta[h, t] for h in model.HarvestCells])
205       
206    # ans3 is T3 (i.e. production costs at origin nodes)
207    ans3 = sum([model.a[h, t] * model.yr[t] * model.A[h] * model.delta[h,t] * model.Q[o, t] for o in model.OriginNodes for h in model.HCellsForOrigin[o]])
208           
209    # ans4 is T4 (i.e. potential road construction cost)
210    ans4 = sum([model.C[i, j, t] * model.gamma[i, j, t] for (i,j) in model.PotentialRoads])
211   
212    # ans5 is T5 (i.e. wood transport cost)
213    ans5 = sum([model.D[i, j, t] * model.f[i, j, t] for (i,j) in model.AllRoads])
214       
215# *****
216# *****
217#    return (model.AnoProfit[t] - (0.9)^(t-1)*(ans1 - ans2 - ans3 - ans4 - ans5)) == 0.0
218# *****
219# *****
220    return (model.AnoProfit[t] - (ans1 - ans2 - ans3 - ans4 - ans5)) == 0.0
221
222model.ComputeStageProfit = Constraint(model.Times, rule=stage_profit_rule)
223
224#
225# Objective
226#
227
228def total_profit_rule(model):
229   return (summation(model.AnoProfit))
230
231model.Production_Profit_Objective = Objective(rule=total_profit_rule, sense=maximize)
Note: See TracBrowser for help on using the repository browser.