1 # Model I Forest Estate Modelling using GLPK/MathProg

2 # Reading and writing dbf files

3

4 # by Noli Sicad --- nsicad@gmail.com

5 # 18 December 2009

6

7 # Forest Management 4th Edition

8 # by Lawrence Davis, K. Norman Johnson, Pete Bettinger, Theodore Howard

9 # Chapter 11 - Daniel Pickett

10 # http://warnell.forestry.uga.edu/Warnell/Bettinger/mgtbook/index.htm

11

12 # Model I Formulation

13

14 /* Note: This is not the full LP model mentioned in the book.

15 Some of the constraints are deliberately omitted in this model for the purpose of clarity.

16

17 The features of MathProg in this example are:

18 * reading and writing dbf from regular dbf files,

19 * reading dbf file (database of shapefile (stands.shp)) (e.g. area parameter),

20 * using the area data in the constraints and

21 * writing dbf file from result of LP model.

22

23 Model I - Harvest Scheduling formulation for Sustainable Forest Management (SFM)

24

25 Features are:

26 * Net Present Value for the objective function (Revenue - Cost)

27 * Harvest Constraints by period - Sustainable Yield per Period

28 * Even-Flow Constraint / Volume - Harvest Flow Constraint - Alpha (1-Apha)

29 * Even-Flow Constraint / Volume - Harvest Flow Constraint - Beta (1 +Beta)

30 * Forest / Land Constraint -- Total Area of the forest

31 * Forest Stand Constraint -- Individual stands

32

33 What is next? -- Forest Mgt Carbon Accounting for Climate Change

34

35 Note: The model file that the data containing in

36 the dbf files is public domain material (so it is compatible with the

37 GNU GPL) and data can be found in

38 http://warnell.forestry.uga.edu/Warnell/Bettinger/mgtbook/index.htm

39

40 # Noli Sicad --- nsicad@gmail.com

41

42 */

43

44 set G_STAND_TYPE; # A, B, C

45

46 set I_CULTURAL_PRES;

47 set J_MGT_YEAR;

48

49 param K_PERIOD;

50 param Forest_Cost{G_STAND_TYPE,I_CULTURAL_PRES, J_MGT_YEAR}; # cost

51

52 param Yield_Table_Vol{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD} >= 0;

53

54

55 param Alpha >= 0;

56 param Beta >= 0;

57

58 param TCost_Table{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD} >= 0;

59

60 param NetRev_Table{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD};

61

62

63 var XForestLand{g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} >= 0;

64

65

66 #reading dbf tables

67 table tab IN "xBASE" "standtype.dbf": G_STAND_TYPE <- [STAND];

68 display G_STAND_TYPE;

69

70

71 table tab2 IN "xBASE" "cultural_pres.dbf": I_CULTURAL_PRES <- [CUL_PRES];

72 display I_CULTURAL_PRES;

73

74 table tab3 IN "xBASE" "mgt_year.dbf": J_MGT_YEAR <- [MGT_YEAR];

75 display J_MGT_YEAR;

76

77 /*

78 param Forest_Cost{G_STAND_TYPE,I_CULTURAL_PRES, J_MGT_YEAR} default 0; # cost

79 */

80

81 set S1, dimen 3;

82 table tab4 IN "xBASE" "Forest_Cost.dbf": S1 <- [STAND, CUL_PRES, MGT_YEAR],Forest_Cost ~FCOST;

83 display Forest_Cost;

84

85 set S2, dimen 4;

86 table tab5 IN "xBASE" "Yield_Table_Vol.dbf": S2 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],Yield_Table_Vol ~YIELD;

87 display Yield_Table_Vol;

88

89 set S3, dimen 4;

90 table tab5 IN "xBASE" "TCost_Table.dbf": S3 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],TCost_Table ~TCOST;

91 display TCost_Table;

92

93

94 set S4, dimen 4;

95 table tab5 IN "xBASE" "NetRev_Table.dbf": S4 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],NetRev_Table ~NETREV;

96 display NetRev_Table;

97

98

99 param MGT;

100

101 param Area_Stand_Indi{g in G_STAND_TYPE, m in 1..MGT} default 0;

102

103 set ST, dimen 2;

104 table tab5 IN "xBASE" "stands.dbf": ST <- [VEG_TYPE, MGT], Area_Stand_Indi ~ACRES;

105 display Area_Stand_Indi;

106

107 param Area_Stand_Type{g in G_STAND_TYPE}:= sum {m in 1..MGT } Area_Stand_Indi[g,m];

108 display Area_Stand_Type;

109

110

111 param Total_Area := sum {g in G_STAND_TYPE, m in 1..MGT } Area_Stand_Indi[g,m];

112 display Total_Area;

113

114 param Harvest_Min_Vol_Period;

115

116

117 var NetPresentValue;

118

119 # Objective function

120 maximize Net_Present_Value: NetPresentValue;

121

122 subject to NPV:

123 NetPresentValue = sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Forest_Cost[g,i,j] * XForestLand[g,i,j];

124

125 # Harvest Constraint by Period

126 subject to Harvest_Period_H {k in 1..K_PERIOD}:

127 sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j] >= Harvest_Min_Vol_Period;

128

129

130 #Even-Flow Constraint / Volume - Harvest Flow Constraint - Alpha

131 subject to Even_Flow_Constaints_Alpha {k in 6..K_PERIOD-1}:

132 (1 - Alpha) * sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j] -

133 sum {g in G_STAND_TYPE,i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k+1] * XForestLand[g,i,j] <= 0;

134

135 # Even-Flow Constraint / Volume - Harvest Flow Constraint - Beta

136 subject to Even_Flow_Constaints_Beta {k in 6..K_PERIOD-1}:

137 (1 + Beta) * sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j] -

138 sum {g in G_STAND_TYPE,i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k+1] * XForestLand[g,i,j] >= 0;

139

140 # Forest / Land Constraints

141 subject to Total_Area_Constraint:

142 sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} XForestLand[g,i,j] <= Total_Area;

143 display Total_Area;

144

145 # Forest / Land Constraints for A B C

146 subject to Area {g in G_STAND_TYPE}:

147 sum {i in I_CULTURAL_PRES,j in J_MGT_YEAR} XForestLand[g,i,j] = Area_Stand_Type[g];

148

149

150

151 solve;

152 #RESULT SECTION

153 printf '#################################\n';

154 printf 'Forest Management Model I - Noli Sicad\n';

155 printf '\n';

156 printf 'Net Present Value = %.2f\n', NetPresentValue;

157 printf '\n';

158

159 printf '\n';

160 printf 'Variables\n';

161 printf 'Stand_Type Age_Class Mgt_Presc Sign Value \n';

162 printf{g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR}:'%5s %10s %11s = %10.2f\n', g,i,j, XForestLand[g,i,j];

163 printf '\n';

164

165 printf 'Constraints\n';

166 printf 'Period Harvest Sign \n';

167 for {k in 1..K_PERIOD} {

168 printf '%5s %10.2f >= %.3f\n', k, sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j], Harvest_Min_Vol_Period;

169 }

170

171 # xbase (dbf) output

172 table Harvest{k in 1..K_PERIOD} OUT "xBASE" "HarvestArea1.dbf" "N(5)N(15,2)" : k ~ Period, (sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j]) ~ H_Area;

173

174 # xbase (dbf) read

175 set S, dimen 2;

176 table tab2 IN "xBASE" "HarvestArea1.dbf": S <- [Period, H_Area];

177 display S;

178

179

180

181

182 printf '\n';

183 printf 'Constraint\n';

184 printf 'Harvest Period\n';

185 printf 'Type AgeClass PrescMgt Period Value\n';

186 printf{g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR, k in 1..K_PERIOD}:'%5s %11s %11s %5s %10.2f\n', g,i,j, k, (Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j]);

187

188

189 printf 'Even_Flow_Constaint_Alpha (1-Alpha)\n';

190 printf 'Period Sign \n';

191 for {k in 6..K_PERIOD-1} {

192 printf "%s %10.2f <= %s\n", k, ((1 - Alpha) * sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k] * XForestLand[g,i,j] - sum {g in G_STAND_TYPE,i in I_CULTURAL_PRES, j in J_MGT_YEAR} Yield_Table_Vol[g,i,j,k+1] * XForestLand[g,i,j]),0;

193 }

194 printf '\n';

195

196

197 # Forest / Land Constraints

198 printf '\n';

199 printf 'Total Area Constraint\n';

200 printf 'Type AgeClass PrescMgt Value Sign Total_Area \n';

201 printf '%5s <= %.3f\n',sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} XForestLand[g,i,j], Total_Area;

202

203 printf 'Area\n';

204 printf 'Area Value Sign Areas_Stand\n';

205 for {g in G_STAND_TYPE} {

206 printf '%5s %10.2f <= %.3f\n', g, sum {i in I_CULTURAL_PRES,j in J_MGT_YEAR} XForestLand[g,i,j], Area_Stand_Type[g];

207 }

208

209

210 #DATA SECTION

211

212 data;

213

214 # Most of the data has been moved to dbf format

215

216 param MGT:=31;

217

218 param K_PERIOD:= 7;

219

220 param Alpha:= 0.20;

221 param Beta:= 0.20;

222

223 param Harvest_Min_Vol_Period:= 12000;

224

225 end;

226