Added infrastructure expansion cost.

GAMS/IntExtOpt.gms

@@ -69,6 +69,7 @@
  PARAMETER timberForestYield(location);
 
  PARAMETER conversionCost(land_cover, land_cover) cost of converting from one land cover to another;
+ PARAMETER infrastructureCost(location);
 
 *$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
 $gdxin %gdxincname%
@@ -78,7 +79,7 @@ $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
 $load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
 $load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
 $load previousLandCoverArea, carbonFluxRate, woodYield, carbonPrice, woodPrice
-$load minimumLandCoverArea, conversionCost, timberForestYield, timberYield
+$load minimumLandCoverArea, conversionCost, timberForestYield, timberYield, infrastructureCost
 $gdxin
 
  SCALAR delta "use to smooth power function see 7.5 www.gams.com/dd/docs/solversconopt.pdf" / 0.00000000001 /;
@@ -117,8 +118,8 @@ $gdxin
  baseCost('pasture') = 0.04;
  otherIntCost('pasture') = 0.8 + otherICost;
 
- SCALAR forestExpansionMaxRate / 1 /;
- SCALAR restrictedDeforestationCost / 10000000 /;
+ SCALAR forestExpansionMaxRate / 0.05 /;
+ SCALAR restrictedDeforestationCost might fail to optimise if set too high / 100 /;
  unrestrictedArea(land_cover, location) = previousLandCoverArea(land_cover, location) - minimumLandCoverArea(land_cover, location);
  woodYield('carbonForest', 'naturalNew', location) = woodYield('carbonForest', 'timberForest', location);
  woodYield('natural', 'naturalNew', location) = woodYield('natural', 'pasture', location);
@@ -188,7 +189,8 @@ $gdxin
        LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
        RESTRICTED_LAND_COVER_CHANGE_CALC(land_cover, location)
        RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location)
-       FOREST_EXPANSION_CONSTRAINT
+       FOREST_EXPANSION_CONSTRAINT(forest, location)
+       FOREST_CONTRACTION_CONSTRAINT(forest, location)
        PREMATURE_DEFORESTATION_COST_CALC(location)
        TOTAL_LAND_COVER_CALC(land_cover, location)
        TOTAL_LAND_COVER_CHANGE_CALC(land_cover, land_cover_after, location)
@@ -273,7 +275,8 @@ $gdxin
 
  NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) .. totalLandCoverChange(exc_natural, 'natural', location) =E= 0;
 
- FOREST_EXPANSION_CONSTRAINT .. sum((non_forest, forest, location), landCoverChange(non_forest, forest, location)) =L= sum((non_forest, location), previousLandCoverArea(non_forest, location)) * forestExpansionMaxRate;
+ FOREST_EXPANSION_CONSTRAINT(forest, location) .. totalLandCoverArea(forest, location) =L= previousLandCoverArea(forest, location) * (1 + forestExpansionMaxRate);
+ FOREST_CONTRACTION_CONSTRAINT(forest, location) .. totalLandCoverArea(forest, location) =G= previousLandCoverArea(forest, location) * (1 - forestExpansionMaxRate);
 
  PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum((forest, land_cover)$[not sameAs(forest, land_cover)], restrictedLandCoverChange(forest, land_cover, location)) * restrictedDeforestationCost;
 
@@ -289,7 +292,7 @@ $gdxin
  CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
 
 * Subtracting cost of not fully grown timberForest reserved under minimumLandCoverArea
- CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * (conversionCost(land_cover_before, land_cover_after) + ((previousLandCoverArea('natural', location) / suitableLandArea(location)) * 0.5 + 0.1)));
+ CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * (conversionCost(land_cover_before, land_cover_after) + infrastructureCost(location)));
 
  COST_EQ .. total_cost =E=
          (

src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java

@@ -116,6 +116,7 @@ public class GamsLocationOptimiser {
 		GAMSParameter seedAndWasteRateP = inDB.addParameter("seedAndWasteRate", 1);	
 		
 		GAMSParameter prevLandCoverP = inDB.addParameter("previousLandCoverArea", 2);
+		GAMSParameter infrastructureCostP = inDB.addParameter("infrastructureCost", 1);
 		
 		double totalAgriLand = 0;
 		double totalSuitable = 0;
@@ -179,6 +180,10 @@ public class GamsLocationOptimiser {
 				v.add(locString);
 				setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(lc), 5);
 			}
+			
+			// Infrastructure cost
+			double infrastructureCost = landUseItem.getLandCoverFract(LandCoverType.NATURAL) * 0.09 + 0.01;
+			setGamsParamValue(infrastructureCostP.addRecord(locString), infrastructureCost, 5);
 				
 		}