source: coopr.pysp/trunk/examples/pysp/forestry/chile18-nb-yr/ReferenceModel.py @ 2037

Last change on this file since 2037 was 2037, checked in by dlwoodr, 10 years ago

add this example from Fernando

  • Property svn:executable set to *
File size: 10.9 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# origin nodes for cell
49model.COriginNodeForCell = Set(model.HarvestCells, within=model.OriginNodes)
50
51# existing roads
52model.ExistingRoads = Set(within=model.Nodes*model.Nodes) 
53
54# potential roads
55model.PotentialRoads = Set(within=model.Nodes*model.Nodes)
56
57# Roads declaring
58#model.AllRoads = model.ExistingRoads & model.PotentialRoads
59model.AllRoads = Set(within=model.Nodes*model.Nodes)
60
61# derived set
62def aplus_arc_set_rule(v, model):
63   return  [(i, j) for (i, j) in model.AllRoads if j == v]
64model.Aplus = Set(model.Nodes, within=model.AllRoads, rule=aplus_arc_set_rule)
65
66# derived set
67def aminus_arc_set_rule(v, model):
68   return  [(i, j) for (i, j) in model.AllRoads if i == v]
69model.Aminus = Set(model.Nodes, within=model.AllRoads, rule=aminus_arc_set_rule)
70
71# derived set
72# no isolated roads can be built
73def connected_potential_roads_rule(k, l, model):
74   return  [ (i, j) for (i, j) in model.PotentialRoads if (i==k and j!=l) or (i==l and j!=k) or (j==k and i!=l) or (j==l and i!=k) ]
75model.ConnectedPotentialRoads = Set( model.Nodes, model.Nodes, within=model.PotentialRoads, rule=connected_potential_roads_rule)
76
77def potential_roads_connected_to_existing_roads_rule(model):
78   set=[]
79   for (k,l) in model.PotentialRoads:
80      if not (k in model.ExitNodes and l in model.OriginNodes) and not (l in model.ExitNodes and k in model.OriginNodes):
81         for x in model.Nodes:
82            if x!=k and x!=l:
83               #print "k,l,x=",k,l,x
84               if (x,k) in model.AllRoads or (x,l) in model.AllRoads or (k,x) in model.AllRoads or (l,x) in model.AllRoads:
85                  #if not (k in model.OriginNodes and (l,x) in model.ExistingRoads) and not (l in model.OriginNodes and (k,x) in model.ExistingRoads):
86                  if ((x,k) in model.ExistingRoads) or ((x,l) in model.ExistingRoads) or ((k,x) in model.ExistingRoads) or ((l,x) in model.ExistingRoads):
87                     if (k,l) not in set:
88                        set.append((k,l))
89                        #print "set=",set
90   return set
91
92model.PotentialRoadsConnectedToExistingRoads = Set(within=model.PotentialRoads, rule=potential_roads_connected_to_existing_roads_rule)
93
94model.DisconnectedPotentialRoads = model.PotentialRoads - model.PotentialRoadsConnectedToExistingRoads
95
96# derived set
97# no isolated lots can be harvest
98def potential_roads_accesing_lot_h(h, model):
99   return  [ (i, j) for (i, j) in model.PotentialRoads if (i in model.COriginNodeForCell[h]) or (j in model.COriginNodeForCell[h]) ]
100model.PotentialRoadsAccesingLotH = Set( model.HarvestCells, within=model.PotentialRoads, rule=potential_roads_accesing_lot_h )
101
102def lots_connected_to_existing_roads(model):
103   set=[]
104   for h in model.HarvestCells:
105      for (i,j) in model.ExistingRoads:
106         if (i in model.COriginNodeForCell[h]) or (j in model.COriginNodeForCell[h]):
107            if h not in set:
108               set.append(h)
109   return set
110
111model.LotsConnectedToExistingRoads = Set(within=model.HarvestCells, rule=lots_connected_to_existing_roads )
112
113model.LotsNotConnectedToExistingRoads =  model.HarvestCells - model.LotsConnectedToExistingRoads
114
115
116#
117# Deterministic Parameters
118#
119
120# productivity of cell h if it is harvested in period t
121model.a = Param(model.HarvestCells, model.Times)
122
123# Area of cell h to be harvested
124model.A = Param(model.HarvestCells)
125
126# flow capacity (smallest big M)
127def Umax_init(model):
128   return (sum([model.a[h, t] * model.A[h] for h in model.HarvestCells for t in model.Times])/2)()
129model.Umax = Param(initialize=Umax_init)
130
131# harvesting cost of one hectare of cell h in time t
132model.P = Param(model.HarvestCells, model.Times)
133
134# unit Production cost at origin o in time t
135model.Q = Param(model.OriginNodes, model.Times)
136
137# construction cost of one road in arc (k,l) in time t
138model.C = Param(model.PotentialRoads, model.Times, default=0.0)
139
140# unit transport cost through arc (k,l) in time t
141model.D = Param(model.AllRoads, model.Times)
142
143#
144# Stochastic Parameters
145#
146
147# sale price at exit s in time period t
148model.R = Param(model.ExitNodes, model.Times)
149
150# upper and lower bounds on wood supplied
151model.Zlb = Param(model.Times)
152model.Zub = Param(model.Times)
153
154# Yield Ratio
155model.yr = Param(model.Times)
156
157#
158# Variables
159#
160
161# If cell h is harvested in period t
162model.delta = Var(model.HarvestCells, model.Times, domain=Boolean)
163
164# If road in arc (k,l) is built in period t
165model.gamma = Var(model.PotentialRoads, model.Times, domain=Boolean)
166
167# Flow of wood thransported through arc (k,l)
168model.f = Var(model.AllRoads, model.Times, domain=NonNegativeReals)
169
170# Supply of wood at exit s in period t
171model.z = Var(model.ExitNodes, model.Times, domain=NonNegativeReals)
172
173# Declared when changing the program for stochastics
174model.AnoProfit = Var(model.Times)
175
176#
177# Constraints
178#
179
180# flow balance at origin nodes
181def origin_flow_bal(o, t, model):
182   ans1 = sum([model.a[h, t] * model.yr[t] * model.A[h] * model.delta[h, t] for h in model.HCellsForOrigin[o]])
183   ans2 = sum([model.f[k, o, t] for (k, o) in model.Aplus[o]])
184   ans3 = sum([model.f[o, k, t] for (o, k) in model.Aminus[o]])
185   return (ans1 + ans2 - ans3) == 0.0
186
187model.EnforceOriginFlowBalance = Constraint(model.OriginNodes, model.Times, rule=origin_flow_bal)
188
189# flow balance at intersection nodes
190def intersection_flow_bal(j, t, model):
191   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
192
193model.EnforceIntersectionFlowBalance = Constraint(model.IntersectionNodes, model.Times, rule=intersection_flow_bal)
194
195# flow balance at destination nodes
196def destination_flow_bal(e, t, model):
197   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
198
199model.EnforceDestinationFlowBalance = Constraint(model.ExitNodes, model.Times, rule=destination_flow_bal)
200
201# explicit divergence or supply equals demand or no inventory between periods
202def divergence_flow_bal(t, model):
203   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
204
205model.EnforceDivergenceFlowBalance = Constraint(model.Times, rule=divergence_flow_bal)
206
207# Wood Production Bounds
208def wood_prod_bounds(t, model):
209   return (model.Zlb[t], sum([model.z[e, t] for e in model.ExitNodes]), model.Zub[t])
210
211model.EnforceWoodProductionBounds = Constraint(model.Times, rule=wood_prod_bounds)
212
213# potential roads flow capacity 1
214def pot_road_flow_capacity(k, l, t, model):
215        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)
216
217model.EnforcePotentialRoadFlowCap = Constraint(model.PotentialRoads, model.Times, rule=pot_road_flow_capacity)
218
219# each potential road can be built no more than once in the time horizon
220def one_potential_road(k, l, model):
221   return (None, sum([model.gamma[k ,l, t] for t in model.Times]), 1.0)
222
223model.EnforceOnePotentialRoad = Constraint(model.PotentialRoads, rule=one_potential_road)
224
225# existing roads flow capacity
226def existing_roads_capacity(k, l, t, model):
227        return (None, model.f[k, l, t], model.Umax)
228
229model.EnforceExistingRoadCapacities = Constraint(model.ExistingRoads, model.Times, rule=existing_roads_capacity)
230
231# additional road bound capacity
232def additional_roads_capacity(k, l, t, model):
233        return ( None, model.f[k, l, t] - sum([model.z[e, t] for e in model.ExitNodes]), 0.0 )
234
235model.EnforceAdditionalRoadCapacities = Constraint(model.AllRoads, model.Times, rule=additional_roads_capacity)
236
237# each cell can be harvested no more than once in the time horizon
238def one_harvest(h, model):
239   return (None, sum([model.delta[h,t] for t in model.Times]), 1.0)
240
241model.EnforceOneHarvest = Constraint(model.HarvestCells, rule=one_harvest)
242
243# no isolated roads can be built
244def road_to_road_trigger(t, k, l, model):
245   return ( None, sum([ model.gamma[ k, l, tt] for tt in model.Times if tt <= t ]) - sum([ model.gamma[ i, j, ttt] for ttt in model.Times if ttt <= t for (i,j) in model.ConnectedPotentialRoads[k,l] ]), 0.0 )
246
247model.EnforceRoadToRoadTrigger = Constraint(model.Times, model.DisconnectedPotentialRoads, rule=road_to_road_trigger)
248
249# no isolated lots can be harvest
250def lot_to_road_trigger(t, h, model):
251   return ( None, sum([ model.delta[ h, tt] for tt in model.Times if tt <= t ]) - sum([ model.gamma[ i, j, ttt] for ttt in model.Times if ttt <= t for (i,j) in model.PotentialRoadsAccesingLotH[h] ]), 0.0 )
252
253model.EnforceLotToRoadTrigger = Constraint(model.Times, model.LotsNotConnectedToExistingRoads, rule=lot_to_road_trigger)
254
255#
256# Stage-specific profit computations
257#
258
259def stage_profit_rule(t, model):
260    # ans1 is T1 (i.e. wood sale revenue)
261    ans1 = sum([model.R[e, t] * model.z[e, t] for e in model.ExitNodes])
262   
263    # ans2 is T2 (i.e. wood harvest cost)
264    ans2 = sum([model.P[h, t] * model.A[h] * model.delta[h, t] for h in model.HarvestCells])
265       
266    # ans3 is T3 (i.e. production costs at origin nodes)
267    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]])
268           
269    # ans4 is T4 (i.e. potential road construction cost)
270    ans4 = sum([model.C[i, j, t] * model.gamma[i, j, t] for (i,j) in model.PotentialRoads])
271   
272    # ans5 is T5 (i.e. wood transport cost)
273    ans5 = sum([model.D[i, j, t] * model.f[i, j, t] for (i,j) in model.AllRoads])
274       
275# *****
276# *****
277#    return (model.AnoProfit[t] - (0.9)^(t-1)*(ans1 - ans2 - ans3 - ans4 - ans5)) == 0.0
278# *****
279# *****
280    return (model.AnoProfit[t] - (ans1 - ans2 - ans3 - ans4 - ans5)) == 0.0
281
282model.ComputeStageProfit = Constraint(model.Times, rule=stage_profit_rule)
283
284#
285# Objective
286#
287
288def total_profit_rule(model):
289   return (summation(model.AnoProfit))
290
291model.Production_Profit_Objective = Objective(rule=total_profit_rule, sense=maximize)
Note: See TracBrowser for help on using the repository browser.