source: coopr.pysp/stable/2.1/examples/pysp/forestry/models-nb-yr/ReferenceModel.py @ 2068

Last change on this file since 2068 was 2068, checked in by wehart, 10 years ago

Merged revisions 1952-2067 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pysp/trunk

........

r1956 | jwatson | 2009-12-02 17:56:53 -0700 (Wed, 02 Dec 2009) | 3 lines


Added --scenario-solver-options and --ef-solver-options options to the "runph" script.

........

r1957 | dlwoodr | 2009-12-03 14:17:35 -0700 (Thu, 03 Dec 2009) | 2 lines


Documentation updates for pysp

........

r1974 | wehart | 2009-12-06 17:20:56 -0700 (Sun, 06 Dec 2009) | 2 lines


Updating PyPI categories

........

r1978 | jwatson | 2009-12-10 21:29:33 -0700 (Thu, 10 Dec 2009) | 3 lines


Eliminated exception-handling logic when loading user-defined extension modules in PH. The range of exceptions is too large, and for debugging purposes, it is more useful to see the raw trace output.

........

r1979 | jwatson | 2009-12-10 22:23:17 -0700 (Thu, 10 Dec 2009) | 5 lines


Biggest enhancement: The efwriter command-line script now has the option to output a CVaR-weighted objective term. Not extensively tested, but behaves sane on a number of small test cases.


Improved exception handling and error diagnostics in both the runph and efwriter scripts.

........

r1985 | jwatson | 2009-12-12 10:45:17 -0700 (Sat, 12 Dec 2009) | 3 lines


Modified PH to only use warm-starts if a solver has the capability!

........

r1998 | jwatson | 2009-12-13 15:17:58 -0700 (Sun, 13 Dec 2009) | 3 lines


Changed references to _component to active_component.

........

r2026 | wehart | 2009-12-21 23:27:06 -0700 (Mon, 21 Dec 2009) | 2 lines


Attempting to update PH. I'm not sure if this works, since I don't know how to test PH.

........

r2029 | jwatson | 2009-12-22 09:52:21 -0700 (Tue, 22 Dec 2009) | 3 lines


Some fixes to the ef writer based on Bill's recent changes to _initialize_constraint.

........

r2035 | jwatson | 2009-12-22 21:10:32 -0700 (Tue, 22 Dec 2009) | 3 lines


Added --scenario-mipgap option to PH script. Added _mipgap attribute to PH object. This attribute is mirrored to the solver plugin at the initiation of each iteration, after any PH extensions have the opportunity to provide a new value to the attribute. It is currently made use of by the WW PH extension.

........

r2037 | dlwoodr | 2009-12-23 14:38:43 -0700 (Wed, 23 Dec 2009) | 2 lines


add this example from Fernando

........

r2038 | dlwoodr | 2009-12-23 14:46:56 -0700 (Wed, 23 Dec 2009) | 3 lines


finish the job: we can now duplicate Fernando's example

........

r2039 | jwatson | 2009-12-23 15:13:54 -0700 (Wed, 23 Dec 2009) | 3 lines


Missed fix with new constraint initialization syntax in PH linearization.

........

r2058 | jwatson | 2009-12-29 10:57:58 -0700 (Tue, 29 Dec 2009) | 3 lines


Missed some _initialize_constraint function calls within the PySP EF writer during the recent switch to the corresponding "add" method.

........

r2059 | jwatson | 2009-12-29 10:58:34 -0700 (Tue, 29 Dec 2009) | 3 lines


Enabling garbage collection by default in PH.

........

r2060 | jwatson | 2009-12-29 10:59:05 -0700 (Tue, 29 Dec 2009) | 3 lines


Elimnating some debug output.

........

r2061 | jwatson | 2009-12-29 11:07:47 -0700 (Tue, 29 Dec 2009) | 3 lines


Fixing some option documentation in PH.

........

r2062 | jwatson | 2009-12-29 12:00:37 -0700 (Tue, 29 Dec 2009) | 3 lines


Added ef-mipgap option to PH scripts.

........

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