examples/dbf/ForestMgt_Model_I_GIS_dbf.mod
changeset 1 c445c931472f
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/examples/dbf/ForestMgt_Model_I_GIS_dbf.mod	Mon Dec 06 13:09:21 2010 +0100
     1.3 @@ -0,0 +1,226 @@
     1.4 +#  Model I Forest Estate Modelling using GLPK/MathProg
     1.5 +#  Reading and writing dbf files
     1.6 +
     1.7 +#  by Noli Sicad --- nsicad@gmail.com
     1.8 +# 18 December 2009
     1.9 +
    1.10 +#  Forest Management 4th Edition 
    1.11 +#  by Lawrence Davis, K. Norman Johnson, Pete Bettinger, Theodore Howard
    1.12 +#  Chapter 11 - Daniel Pickett 
    1.13 +#  http://warnell.forestry.uga.edu/Warnell/Bettinger/mgtbook/index.htm
    1.14 +
    1.15 +#  Model I Formulation
    1.16 +
    1.17 +/*  Note: This is not the full LP model mentioned in the book.
    1.18 +Some of the constraints are deliberately omitted in this model for the purpose of clarity.
    1.19 +
    1.20 +The features of MathProg in this example are:
    1.21 +* reading and writing dbf from regular dbf files,
    1.22 +* reading dbf file (database of shapefile (stands.shp)) (e.g. area parameter),
    1.23 +* using the area data in the constraints and
    1.24 +* writing dbf file from result of LP model.
    1.25 +
    1.26 +Model I - Harvest Scheduling formulation for Sustainable Forest Management (SFM)
    1.27 +
    1.28 +Features are:
    1.29 +* Net Present Value for the objective function (Revenue - Cost)
    1.30 +* Harvest Constraints by period - Sustainable Yield per Period
    1.31 +* Even-Flow Constraint / Volume - Harvest Flow Constraint -  Alpha (1-Apha)
    1.32 +* Even-Flow Constraint / Volume - Harvest Flow Constraint - Beta  (1 +Beta)
    1.33 +* Forest / Land Constraint -- Total Area of the forest
    1.34 +* Forest Stand Constraint  -- Individual stands
    1.35 +
    1.36 +What is next? -- Forest Mgt Carbon Accounting for Climate Change
    1.37 +
    1.38 +Note: The model file that the data containing in
    1.39 +the dbf files is public domain material (so it is compatible with the
    1.40 +GNU GPL) and data can be found in 
    1.41 +http://warnell.forestry.uga.edu/Warnell/Bettinger/mgtbook/index.htm
    1.42 +
    1.43 +# Noli Sicad --- nsicad@gmail.com
    1.44 +
    1.45 +*/
    1.46 +
    1.47 +set G_STAND_TYPE; # A, B, C
    1.48 +
    1.49 +set I_CULTURAL_PRES; 
    1.50 +set J_MGT_YEAR; 
    1.51 +
    1.52 +param K_PERIOD;
    1.53 +param Forest_Cost{G_STAND_TYPE,I_CULTURAL_PRES, J_MGT_YEAR}; # cost
    1.54 +
    1.55 +param Yield_Table_Vol{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD} >= 0;
    1.56 +
    1.57 +
    1.58 +param Alpha >= 0;
    1.59 +param Beta >= 0;
    1.60 +
    1.61 +param TCost_Table{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD} >= 0;
    1.62 +
    1.63 +param NetRev_Table{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD};
    1.64 +
    1.65 +
    1.66 +var XForestLand{g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} >= 0;
    1.67 +
    1.68 +
    1.69 +#reading dbf tables
    1.70 +table tab IN "xBASE" "standtype.dbf": G_STAND_TYPE <- [STAND];
    1.71 +display G_STAND_TYPE;
    1.72 +
    1.73 +
    1.74 +table tab2 IN "xBASE" "cultural_pres.dbf": I_CULTURAL_PRES <- [CUL_PRES];
    1.75 +display I_CULTURAL_PRES;
    1.76 +
    1.77 +table tab3 IN "xBASE" "mgt_year.dbf": J_MGT_YEAR <- [MGT_YEAR];
    1.78 +display J_MGT_YEAR;
    1.79 +
    1.80 +/*
    1.81 +param Forest_Cost{G_STAND_TYPE,I_CULTURAL_PRES, J_MGT_YEAR} default 0; # cost
    1.82 +*/
    1.83 +
    1.84 +set S1, dimen 3;
    1.85 +table tab4 IN "xBASE" "Forest_Cost.dbf": S1 <- [STAND, CUL_PRES, MGT_YEAR],Forest_Cost ~FCOST;
    1.86 +display Forest_Cost;
    1.87 +
    1.88 +set S2, dimen 4;
    1.89 +table tab5 IN "xBASE" "Yield_Table_Vol.dbf": S2 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],Yield_Table_Vol ~YIELD;
    1.90 +display Yield_Table_Vol;
    1.91 +
    1.92 +set S3, dimen 4;
    1.93 +table tab5 IN "xBASE" "TCost_Table.dbf": S3 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],TCost_Table ~TCOST;
    1.94 +display TCost_Table;
    1.95 +
    1.96 +
    1.97 +set S4, dimen 4;
    1.98 +table tab5 IN "xBASE" "NetRev_Table.dbf": S4 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],NetRev_Table ~NETREV;
    1.99 +display NetRev_Table;
   1.100 +
   1.101 +
   1.102 +param MGT;
   1.103 +
   1.104 +param Area_Stand_Indi{g in G_STAND_TYPE, m in 1..MGT} default 0; 
   1.105 +
   1.106 +set ST, dimen 2;
   1.107 +table tab5 IN "xBASE" "stands.dbf": ST <- [VEG_TYPE, MGT], Area_Stand_Indi ~ACRES;
   1.108 +display Area_Stand_Indi;
   1.109 +
   1.110 +param Area_Stand_Type{g in G_STAND_TYPE}:= sum {m in 1..MGT } Area_Stand_Indi[g,m];
   1.111 +display Area_Stand_Type;
   1.112 +
   1.113 +
   1.114 +param Total_Area := sum {g in G_STAND_TYPE, m in 1..MGT } Area_Stand_Indi[g,m];
   1.115 +display Total_Area;
   1.116 +
   1.117 +param Harvest_Min_Vol_Period;
   1.118 +
   1.119 +
   1.120 +var NetPresentValue;
   1.121 +
   1.122 +# Objective function
   1.123 +maximize Net_Present_Value: NetPresentValue;
   1.124 +
   1.125 +subject to NPV:
   1.126 +   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];
   1.127 +
   1.128 +# Harvest Constraint by Period
   1.129 +subject to Harvest_Period_H {k in 1..K_PERIOD}:
   1.130 +   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;
   1.131 +   
   1.132 +
   1.133 +#Even-Flow Constraint / Volume - Harvest Flow Constraint - Alpha
   1.134 +subject to Even_Flow_Constaints_Alpha {k in 6..K_PERIOD-1}:
   1.135 +    (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] -
   1.136 +    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;
   1.137 +
   1.138 +# Even-Flow Constraint / Volume - Harvest Flow Constraint - Beta
   1.139 +subject to Even_Flow_Constaints_Beta {k in 6..K_PERIOD-1}:
   1.140 +    (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] -
   1.141 +    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;
   1.142 +   
   1.143 +# Forest / Land Constraints
   1.144 +subject to Total_Area_Constraint: 
   1.145 +  sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} XForestLand[g,i,j] <= Total_Area;
   1.146 +display Total_Area;   
   1.147 +
   1.148 +# Forest / Land Constraints for A B C
   1.149 +subject to Area {g in G_STAND_TYPE}:
   1.150 +   sum {i in I_CULTURAL_PRES,j in J_MGT_YEAR} XForestLand[g,i,j] = Area_Stand_Type[g];
   1.151 +
   1.152 +
   1.153 +
   1.154 +solve;
   1.155 +#RESULT SECTION
   1.156 +printf '#################################\n';
   1.157 +printf 'Forest Management Model I - Noli Sicad\n';
   1.158 +printf '\n';
   1.159 +printf 'Net Present Value = %.2f\n', NetPresentValue;
   1.160 +printf '\n';
   1.161 +
   1.162 +printf '\n';
   1.163 +printf 'Variables\n';
   1.164 +printf 'Stand_Type  Age_Class  Mgt_Presc  Sign Value \n';
   1.165 +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]; 
   1.166 +printf '\n';
   1.167 +
   1.168 +printf 'Constraints\n';
   1.169 +printf 'Period Harvest Sign \n';
   1.170 +for {k in 1..K_PERIOD} {
   1.171 + 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;
   1.172 +   }
   1.173 +
   1.174 +# xbase (dbf) output
   1.175 +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;
   1.176 +
   1.177 +# xbase (dbf) read
   1.178 +set S, dimen 2;
   1.179 +table tab2 IN "xBASE" "HarvestArea1.dbf": S <- [Period, H_Area];
   1.180 +display S;
   1.181 +
   1.182 +
   1.183 +
   1.184 +
   1.185 +printf '\n';
   1.186 +printf 'Constraint\n';
   1.187 +printf 'Harvest Period\n';
   1.188 +printf 'Type AgeClass  PrescMgt Period    Value\n';
   1.189 +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]); 
   1.190 +
   1.191 +
   1.192 +printf 'Even_Flow_Constaint_Alpha (1-Alpha)\n';
   1.193 +printf 'Period Sign \n';
   1.194 +for {k in 6..K_PERIOD-1} {
   1.195 +   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;
   1.196 +  }
   1.197 +printf '\n';
   1.198 +
   1.199 +
   1.200 +# Forest / Land Constraints
   1.201 +printf '\n';  
   1.202 +printf 'Total Area Constraint\n';
   1.203 +printf 'Type AgeClass  PrescMgt  Value Sign Total_Area \n';
   1.204 +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;
   1.205 +
   1.206 +printf 'Area\n';
   1.207 +printf 'Area Value Sign Areas_Stand\n';
   1.208 +for {g in G_STAND_TYPE} {
   1.209 +  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];
   1.210 +   }
   1.211 +
   1.212 +
   1.213 +#DATA SECTION 
   1.214 +      
   1.215 +data;
   1.216 +
   1.217 +# Most of the data has been moved to dbf format
   1.218 +
   1.219 +param MGT:=31;
   1.220 +
   1.221 +param K_PERIOD:= 7;
   1.222 +
   1.223 +param Alpha:= 0.20;
   1.224 +param Beta:= 0.20;
   1.225 +
   1.226 +param Harvest_Min_Vol_Period:= 12000;
   1.227 +
   1.228 +end;
   1.229 +