alpar@1: #  Model I Forest Estate Modelling using GLPK/MathProg
alpar@1: #  Reading and writing dbf files
alpar@1: 
alpar@1: #  by Noli Sicad --- nsicad@gmail.com
alpar@1: # 18 December 2009
alpar@1: 
alpar@1: #  Forest Management 4th Edition 
alpar@1: #  by Lawrence Davis, K. Norman Johnson, Pete Bettinger, Theodore Howard
alpar@1: #  Chapter 11 - Daniel Pickett 
alpar@1: #  http://warnell.forestry.uga.edu/Warnell/Bettinger/mgtbook/index.htm
alpar@1: 
alpar@1: #  Model I Formulation
alpar@1: 
alpar@1: /*  Note: This is not the full LP model mentioned in the book.
alpar@1: Some of the constraints are deliberately omitted in this model for the purpose of clarity.
alpar@1: 
alpar@1: The features of MathProg in this example are:
alpar@1: * reading and writing dbf from regular dbf files,
alpar@1: * reading dbf file (database of shapefile (stands.shp)) (e.g. area parameter),
alpar@1: * using the area data in the constraints and
alpar@1: * writing dbf file from result of LP model.
alpar@1: 
alpar@1: Model I - Harvest Scheduling formulation for Sustainable Forest Management (SFM)
alpar@1: 
alpar@1: Features are:
alpar@1: * Net Present Value for the objective function (Revenue - Cost)
alpar@1: * Harvest Constraints by period - Sustainable Yield per Period
alpar@1: * Even-Flow Constraint / Volume - Harvest Flow Constraint -  Alpha (1-Apha)
alpar@1: * Even-Flow Constraint / Volume - Harvest Flow Constraint - Beta  (1 +Beta)
alpar@1: * Forest / Land Constraint -- Total Area of the forest
alpar@1: * Forest Stand Constraint  -- Individual stands
alpar@1: 
alpar@1: What is next? -- Forest Mgt Carbon Accounting for Climate Change
alpar@1: 
alpar@1: Note: The model file that the data containing in
alpar@1: the dbf files is public domain material (so it is compatible with the
alpar@1: GNU GPL) and data can be found in 
alpar@1: http://warnell.forestry.uga.edu/Warnell/Bettinger/mgtbook/index.htm
alpar@1: 
alpar@1: # Noli Sicad --- nsicad@gmail.com
alpar@1: 
alpar@1: */
alpar@1: 
alpar@1: set G_STAND_TYPE; # A, B, C
alpar@1: 
alpar@1: set I_CULTURAL_PRES; 
alpar@1: set J_MGT_YEAR; 
alpar@1: 
alpar@1: param K_PERIOD;
alpar@1: param Forest_Cost{G_STAND_TYPE,I_CULTURAL_PRES, J_MGT_YEAR}; # cost
alpar@1: 
alpar@1: param Yield_Table_Vol{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD} >= 0;
alpar@1: 
alpar@1: 
alpar@1: param Alpha >= 0;
alpar@1: param Beta >= 0;
alpar@1: 
alpar@1: param TCost_Table{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD} >= 0;
alpar@1: 
alpar@1: param NetRev_Table{G_STAND_TYPE, I_CULTURAL_PRES, J_MGT_YEAR, 1..K_PERIOD};
alpar@1: 
alpar@1: 
alpar@1: var XForestLand{g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} >= 0;
alpar@1: 
alpar@1: 
alpar@1: #reading dbf tables
alpar@1: table tab IN "xBASE" "standtype.dbf": G_STAND_TYPE <- [STAND];
alpar@1: display G_STAND_TYPE;
alpar@1: 
alpar@1: 
alpar@1: table tab2 IN "xBASE" "cultural_pres.dbf": I_CULTURAL_PRES <- [CUL_PRES];
alpar@1: display I_CULTURAL_PRES;
alpar@1: 
alpar@1: table tab3 IN "xBASE" "mgt_year.dbf": J_MGT_YEAR <- [MGT_YEAR];
alpar@1: display J_MGT_YEAR;
alpar@1: 
alpar@1: /*
alpar@1: param Forest_Cost{G_STAND_TYPE,I_CULTURAL_PRES, J_MGT_YEAR} default 0; # cost
alpar@1: */
alpar@1: 
alpar@1: set S1, dimen 3;
alpar@1: table tab4 IN "xBASE" "Forest_Cost.dbf": S1 <- [STAND, CUL_PRES, MGT_YEAR],Forest_Cost ~FCOST;
alpar@1: display Forest_Cost;
alpar@1: 
alpar@1: set S2, dimen 4;
alpar@1: table tab5 IN "xBASE" "Yield_Table_Vol.dbf": S2 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],Yield_Table_Vol ~YIELD;
alpar@1: display Yield_Table_Vol;
alpar@1: 
alpar@1: set S3, dimen 4;
alpar@1: table tab5 IN "xBASE" "TCost_Table.dbf": S3 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],TCost_Table ~TCOST;
alpar@1: display TCost_Table;
alpar@1: 
alpar@1: 
alpar@1: set S4, dimen 4;
alpar@1: table tab5 IN "xBASE" "NetRev_Table.dbf": S4 <- [STAND, CUL_PRES, MGT_YEAR, PERIOD],NetRev_Table ~NETREV;
alpar@1: display NetRev_Table;
alpar@1: 
alpar@1: 
alpar@1: param MGT;
alpar@1: 
alpar@1: param Area_Stand_Indi{g in G_STAND_TYPE, m in 1..MGT} default 0; 
alpar@1: 
alpar@1: set ST, dimen 2;
alpar@1: table tab5 IN "xBASE" "stands.dbf": ST <- [VEG_TYPE, MGT], Area_Stand_Indi ~ACRES;
alpar@1: display Area_Stand_Indi;
alpar@1: 
alpar@1: param Area_Stand_Type{g in G_STAND_TYPE}:= sum {m in 1..MGT } Area_Stand_Indi[g,m];
alpar@1: display Area_Stand_Type;
alpar@1: 
alpar@1: 
alpar@1: param Total_Area := sum {g in G_STAND_TYPE, m in 1..MGT } Area_Stand_Indi[g,m];
alpar@1: display Total_Area;
alpar@1: 
alpar@1: param Harvest_Min_Vol_Period;
alpar@1: 
alpar@1: 
alpar@1: var NetPresentValue;
alpar@1: 
alpar@1: # Objective function
alpar@1: maximize Net_Present_Value: NetPresentValue;
alpar@1: 
alpar@1: subject to NPV:
alpar@1:    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];
alpar@1: 
alpar@1: # Harvest Constraint by Period
alpar@1: subject to Harvest_Period_H {k in 1..K_PERIOD}:
alpar@1:    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;
alpar@1:    
alpar@1: 
alpar@1: #Even-Flow Constraint / Volume - Harvest Flow Constraint - Alpha
alpar@1: subject to Even_Flow_Constaints_Alpha {k in 6..K_PERIOD-1}:
alpar@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] -
alpar@1:     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;
alpar@1: 
alpar@1: # Even-Flow Constraint / Volume - Harvest Flow Constraint - Beta
alpar@1: subject to Even_Flow_Constaints_Beta {k in 6..K_PERIOD-1}:
alpar@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] -
alpar@1:     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;
alpar@1:    
alpar@1: # Forest / Land Constraints
alpar@1: subject to Total_Area_Constraint: 
alpar@1:   sum {g in G_STAND_TYPE, i in I_CULTURAL_PRES, j in J_MGT_YEAR} XForestLand[g,i,j] <= Total_Area;
alpar@1: display Total_Area;   
alpar@1: 
alpar@1: # Forest / Land Constraints for A B C
alpar@1: subject to Area {g in G_STAND_TYPE}:
alpar@1:    sum {i in I_CULTURAL_PRES,j in J_MGT_YEAR} XForestLand[g,i,j] = Area_Stand_Type[g];
alpar@1: 
alpar@1: 
alpar@1: 
alpar@1: solve;
alpar@1: #RESULT SECTION
alpar@1: printf '#################################\n';
alpar@1: printf 'Forest Management Model I - Noli Sicad\n';
alpar@1: printf '\n';
alpar@1: printf 'Net Present Value = %.2f\n', NetPresentValue;
alpar@1: printf '\n';
alpar@1: 
alpar@1: printf '\n';
alpar@1: printf 'Variables\n';
alpar@1: printf 'Stand_Type  Age_Class  Mgt_Presc  Sign Value \n';
alpar@1: 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]; 
alpar@1: printf '\n';
alpar@1: 
alpar@1: printf 'Constraints\n';
alpar@1: printf 'Period Harvest Sign \n';
alpar@1: for {k in 1..K_PERIOD} {
alpar@1:  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;
alpar@1:    }
alpar@1: 
alpar@1: # xbase (dbf) output
alpar@1: 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;
alpar@1: 
alpar@1: # xbase (dbf) read
alpar@1: set S, dimen 2;
alpar@1: table tab2 IN "xBASE" "HarvestArea1.dbf": S <- [Period, H_Area];
alpar@1: display S;
alpar@1: 
alpar@1: 
alpar@1: 
alpar@1: 
alpar@1: printf '\n';
alpar@1: printf 'Constraint\n';
alpar@1: printf 'Harvest Period\n';
alpar@1: printf 'Type AgeClass  PrescMgt Period    Value\n';
alpar@1: 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]); 
alpar@1: 
alpar@1: 
alpar@1: printf 'Even_Flow_Constaint_Alpha (1-Alpha)\n';
alpar@1: printf 'Period Sign \n';
alpar@1: for {k in 6..K_PERIOD-1} {
alpar@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;
alpar@1:   }
alpar@1: printf '\n';
alpar@1: 
alpar@1: 
alpar@1: # Forest / Land Constraints
alpar@1: printf '\n';  
alpar@1: printf 'Total Area Constraint\n';
alpar@1: printf 'Type AgeClass  PrescMgt  Value Sign Total_Area \n';
alpar@1: 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;
alpar@1: 
alpar@1: printf 'Area\n';
alpar@1: printf 'Area Value Sign Areas_Stand\n';
alpar@1: for {g in G_STAND_TYPE} {
alpar@1:   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];
alpar@1:    }
alpar@1: 
alpar@1: 
alpar@1: #DATA SECTION 
alpar@1:       
alpar@1: data;
alpar@1: 
alpar@1: # Most of the data has been moved to dbf format
alpar@1: 
alpar@1: param MGT:=31;
alpar@1: 
alpar@1: param K_PERIOD:= 7;
alpar@1: 
alpar@1: param Alpha:= 0.20;
alpar@1: param Beta:= 0.20;
alpar@1: 
alpar@1: param Harvest_Min_Vol_Period:= 12000;
alpar@1: 
alpar@1: end;
alpar@1: