#  Model I Forest Estate Modelling using GLPK/MathProg

 #  Reading and writing dbf files

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

 # 18 December 2009

 #  Forest Management 4th Edition 

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

 #  Chapter 11 - Daniel Pickett 

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

 #  Model I Formulation

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

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

 The features of MathProg in this example are:

 * reading and writing dbf from regular dbf files,

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

 * using the area data in the constraints and

 * writing dbf file from result of LP model.

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

 Features are:

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

 * Harvest Constraints by period - Sustainable Yield per Period

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

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

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

 * Forest Stand Constraint  -- Individual stands

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

 Note: The model file that the data containing in

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

 GNU GPL) and data can be found in 

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

 # Noli Sicad --- nsicad@gmail.com

 */

 set G_STAND_TYPE; # A, B, C

 set I_CULTURAL_PRES; 

 set J_MGT_YEAR; 

 param K_PERIOD;

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

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

 param Alpha >= 0;

 param Beta >= 0;

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

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

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

 #reading dbf tables

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

 display G_STAND_TYPE;

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

 display I_CULTURAL_PRES;

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

 display J_MGT_YEAR;

 /*

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

 */

 set S1, dimen 3;

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

 display Forest_Cost;

 set S2, dimen 4;

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

 display Yield_Table_Vol;

 set S3, dimen 4;

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

 display TCost_Table;

 set S4, dimen 4;

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

 display NetRev_Table;

 param MGT;

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

 set ST, dimen 2;

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

 display Area_Stand_Indi;

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

 display Area_Stand_Type;

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

 display Total_Area;

 param Harvest_Min_Vol_Period;

 var NetPresentValue;

 # Objective function

 maximize Net_Present_Value: NetPresentValue;

 subject to NPV:

    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];

 # Harvest Constraint by Period

 subject to Harvest_Period_H {k in 1..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] >= Harvest_Min_Vol_Period;

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

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

     (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;

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

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

     (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] -

     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;

 # Forest / Land Constraints

 subject to Total_Area_Constraint: 

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

 display Total_Area;   

 # Forest / Land Constraints for A B C

 subject to Area {g in G_STAND_TYPE}:

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

 solve;

 #RESULT SECTION

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

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

 printf '\n';

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

 printf '\n';

 printf '\n';

 printf 'Variables\n';

 printf 'Stand_Type  Age_Class  Mgt_Presc  Sign Value \n';

 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]; 

 printf '\n';

 printf 'Constraints\n';

 printf 'Period Harvest Sign \n';

 for {k in 1..K_PERIOD} {

  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;

    }

 # xbase (dbf) output

 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;

 # xbase (dbf) read

 set S, dimen 2;

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

 display S;

 printf '\n';

 printf 'Constraint\n';

 printf 'Harvest Period\n';

 printf 'Type AgeClass  PrescMgt Period    Value\n';

 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]); 

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

 printf 'Period Sign \n';

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

    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;

   }

 printf '\n';

 # Forest / Land Constraints

 printf '\n';  

 printf 'Total Area Constraint\n';

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

 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;

 printf 'Area\n';

 printf 'Area Value Sign Areas_Stand\n';

 for {g in G_STAND_TYPE} {

   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];

    }

 #DATA SECTION 

 data;

 # Most of the data has been moved to dbf format

 param MGT:=31;

 param K_PERIOD:= 7;

 param Alpha:= 0.20;

 param Beta:= 0.20;

 param Harvest_Min_Vol_Period:= 12000;

 end;