Refactored GAMS code.

GAMS/IntExtOpt.gms

@@ -58,6 +58,7 @@
 
  PARAMETER previousLandCoverArea(land_cover, location) land cover area in Mha;
  PARAMETER minimumLandCoverArea(land_cover, location) minimum land cover area to constrain conversion in Mha;
+ PARAMETER unrestrictedArea(land_cover, location)
  PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location) carbon flux - MtC-eq per Mha;
  PARAMETER woodYield(land_cover_before, land_cover_after, location)      wood yield - MtC-eq per Mha;
  SCALAR carbonPrice                                 price of carbon - $1000 per tonne;
@@ -115,6 +116,7 @@ $gdxin
 
  SCALAR forestExpansionMaxRate / 1 /;
  SCALAR restrictedDeforestationCost / 10000000 /;
+ unrestrictedArea(land_cover, location) = previousLandCoverArea(land_cover, location) - minimumLandCoverArea(land_cover, location);
 
  VARIABLES
        cropArea(crop, location)             total area for crops - Mha
@@ -130,26 +132,27 @@ $gdxin
        net_supply(crop)                   supply after exports and feed
        landCoverArea(land_cover, location)  land cover area in Mha
        landCoverChange(land_cover_before, land_cover_after, location) land cover change
-       deforestedArea(forest, location) area of forest cut down
-       reservedAreaDeforested(forest, location) area of reserved (minimumLandCoverArea) forest cut down
+       restrictedLandCoverArea(land_cover, location)
+       restrictedLandCoverChange(land_cover, land_cover_after, location)
        prematureDeforestationCost
+       totalLandCoverArea(land_cover, location)
+       totalLandCoverChange(land_cover, land_cover_after, location)
        woodHarvest(location)              total wood harvested
        carbonFlux(location)               total carbon flux
        timberHarvest(location)
        timberForwardContract(location)
-*       potentialWoodYield(location)       potential wood yield for timber forests planted
        totalConversionCost(location)
        totalFeedDM
-*       A                                  "artificial variable https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
+       A                                  "artificial variable for debugging https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
        total_cost                         total cost of domestic supply including net imports;
 
-
  POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
                    totalFeedDM,
-                   deforestedArea, reservedAreaDeforested, prematureDeforestationCost,
-                   landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost;
+                   prematureDeforestationCost,
+                   landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost,
+                   restrictedLandCoverArea, restrictedLandCoverChange, totalLandCoverArea, totalLandCoverChange;
 
-* POSITIVE VARIABLE A;
+ POSITIVE VARIABLE A;
 
  EQUATIONS
        UNIT_COST_EQ(crop, location)                     cost per area - $1000 per ha or $billion per Mha
@@ -165,7 +168,7 @@ $gdxin
        MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location)   constraint on maximum irrigation intensity
        MAX_OTHER_INTENSITY_CONSTRAINT(crop, location)
        SETASIDE_AREA_CALC(location)
-       TOTAL_LAND_CHANGE_CONSTRAINT(location)           constraint on suitable land use
+*       TOTAL_LAND_CHANGE_CONSTRAINT(location)           constraint on suitable land use
        MAX_NET_IMPORT_CONSTRAINT(import_crop)           constraint on max net imports
        MIN_NET_IMPORT_CONSTRAINT(import_crop)           constraint on min net imports
        PASTURE_IMPORT_CONSTRAINT                        constraint to not import pasture
@@ -174,21 +177,19 @@ $gdxin
        NET_SUPPLY_EQ(crop)                              calc net supply for crops
        CROPLAND_LAND_COVER_CALC(location) cropland area equals sum of all crop areas
        PASTURE_LAND_COVER_CALC(location)  pasture area (land cover) equals pasture area (land use)
-*       MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) constraint on land cover conversion
        LAND_COVER_CHANGE_CALC(land_cover, location)
        LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
-*       LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) calculate no change area (e.g. cropland to cropland)
+       RESTRICTED_LAND_COVER_CHANGE_CALC(land_cover, location)
+       RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location)
        FOREST_EXPANSION_CONSTRAINT
-       DEFORESTED_AREA_CALC(forest, location)
-       RESERVED_AREA_DEFORESTED_CALC(forest, location)
        PREMATURE_DEFORESTATION_COST_CALC
-*       UNMANAGED_FOREST_CONSTRAINT(location)                unmanagedForest cannot be created
+       TOTAL_LAND_COVER_CALC(land_cover, location)
+       TOTAL_LAND_COVER_CHANGE_CALC(land_cover, land_cover_after, location)
+*       TOTAL_LAND_COVER_CONSTRAINT(location)
        WOOD_HARVEST_CALC(location)                          calc wood harvested
        CARBON_FLUX_CALC(location)                           calc carbon flux
        TIMBER_HARVEST_CALC(location)
        TIMBER_FORWARD_CONTRACT_CALC(location)
-*       POTENTIAL_WOOD_YIELD_CALC(location)
-*       TIMBER_FOREST_CHANGE_CONSTRAINT(location)                   forcing mature timber forest to be cut down
        CONVERSION_COST(location)
        COST_EQ                                        total cost objective function;
 
@@ -227,10 +228,12 @@ $gdxin
 
  SETASIDE_AREA_CALC(location) .. cropArea('setaside', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) * setAsideRate;
 
-* TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate) + cropArea('pasture', location);
+ TOTAL_LAND_COVER_CALC(land_cover, location) .. totalLandCoverArea(land_cover, location) =E= landCoverArea(land_cover, location) + restrictedLandCoverArea(land_cover, location);
+
+* TOTAL_LAND_COVER_CONSTRAINT(location) .. sum(land_cover, totalLandCoverArea(land_cover, location)) =E= sum(land_cover, previousLandCoverArea(land_cover, location));
 
-* Needs inequality due to floating point errors
- TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(land_cover, landCoverArea(land_cover, location));
+* Needs inequality due to floating point errors but in theory should be equal
+* TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(land_cover, totalLandCoverArea(land_cover, location));
 
  MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
  MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
@@ -239,59 +242,45 @@ $gdxin
 
  IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) * (1.0 - unhandledCropRate) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * cropArea(crop, location));
 
- CROPLAND_LAND_COVER_CALC(location) .. landCoverArea('cropland', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate);
+ CROPLAND_LAND_COVER_CALC(location) .. totalLandCoverArea('cropland', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate);
 
- PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
+ PASTURE_LAND_COVER_CALC(location) .. totalLandCoverArea('pasture', location) =E= cropArea('pasture', location);
 
 
- LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= previousLandCoverArea(land_cover, location) +
+ LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= unrestrictedArea(land_cover, location) +
                                                                                          sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) -
                                                                                          sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
 
-* landCoverChange(A, A, location) defined as unchanged land cover
- LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= previousLandCoverArea(land_cover, location);
-
-* is this needed or can it be handled by keeping track of wood yields?
-* TIMBER_FOREST_CHANGE_CONSTRAINT(location) .. sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location)) =G= matureTimberForestArea(location);
-
- FOREST_EXPANSION_CONSTRAINT .. sum((non_forest, forest, location), landCoverChange(non_forest, forest, location)) =L= sum((non_forest, location), previousLandCoverArea(non_forest, location)) * forestExpansionMaxRate;
+ RESTRICTED_LAND_COVER_CHANGE_CALC(land_cover, location) .. restrictedLandCoverArea(land_cover, location) =E= minimumLandCoverArea(land_cover, location) +
+                                                                                         sum(land_cover_before, restrictedLandCoverChange(land_cover_before, land_cover, location)) -
+                                                                                         sum(land_cover_after, restrictedLandCoverChange(land_cover, land_cover_after, location));
 
-* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =G= minimumLandCoverArea(land_cover, location);
+ TOTAL_LAND_COVER_CHANGE_CALC(land_cover, land_cover_after, location) .. totalLandCoverChange(land_cover, land_cover_after, location) =E= landCoverChange(land_cover, land_cover_after, location) +
+                                                                                                                                          restrictedLandCoverChange(land_cover, land_cover_after, location);
 
-* PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum(forest, (sum((land_cover_after)$[not sameAs(forest, land_cover_after)], landCoverChange(forest, land_cover_after, location)) -
-*                                                         previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)));
-
- DEFORESTED_AREA_CALC(forest, location) .. deforestedArea(forest, location) =E= sum((land_cover_after)$[not sameAs(forest, land_cover_after)], landCoverChange(forest, land_cover_after, location));
-
-* Equivalent to MAX(0, deforestedArea(forest, location) - minimumLandCoverArea(forest, location)) https://www.gams.com/latest/docs/UG_NLP_GoodFormulations.html#UG_NLP_GoodFormulations_ReformulatingDNLPModels
- RESERVED_AREA_DEFORESTED_CALC(forest, location) .. reservedAreaDeforested(forest, location) =E= [(deforestedArea(forest, location) - [previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)] + delta) +
-                                                                                          sqrt(sqr(deforestedArea(forest, location) - [previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)] + delta)) +
-                                                                                               sqr(delta)] / 2;
+* landCoverChange(A, A, location) defined as unchanged land cover
+ LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= unrestrictedArea(land_cover, location);
 
- PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, location), reservedAreaDeforested(forest, location)) * restrictedDeforestationCost;
+ RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, restrictedLandCoverChange(land_cover, land_cover_after, location)) =E= minimumLandCoverArea(land_cover, location);
 
-* UNMANAGED_FOREST_CONSTRAINT(location) .. sum((land_cover)$[not sameAs(land_cover, 'unmanagedForest')], landCoverChange(land_cover, 'unmanagedForest', 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;
 
-* POTENTIAL_WOOD_YIELD_CALC(location) .. potentialWoodYield(location) =E= sum(land_cover_before, landCoverChange(land_cover_before, 'timberForest', location) * woodYield(land_cover_before, 'timberForest', location)) /
-*                                                                         sum(land_cover_before, landCoverChange(land_cover_before, 'timberForest', location));
+ PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, land_cover, location)$[not sameAs(forest, land_cover)], restrictedLandCoverChange(forest, land_cover, location)) * restrictedDeforestationCost;
 
 * Timber from land cover change
  WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((non_timber_forest_before, non_timber_forest_after), landCoverChange(non_timber_forest_before, non_timber_forest_after, location) * woodYield(non_timber_forest_before, non_timber_forest_after, location));
 
 * Timber from timberForest rotation
- TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E=  (sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location)) - minimumLandCoverArea('timberForest', location)) * timberForestYield(location);
+ TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E=  (sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location))) * timberForestYield(location);
 
 * Future timber from newly planted timberForest
  TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(non_timber_forest_before, landCoverChange(non_timber_forest_before, 'timberForest', location) * woodYield(non_timber_forest_before, 'timberForest', location)) +
-                                                                               (landCoverChange('timberForest', 'timberForest', location) - minimumLandCoverArea('timberForest', location)) * woodYield('timberForest', 'timberForest', location);
+                                                                               (landCoverChange('timberForest', 'timberForest', location)) * woodYield('timberForest', 'timberForest', location);
 
- CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
+ 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), landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after));
-
-* CROPLAND_INCREASE_CONSTRAINT(location) .. cropIncrease(location)/(1.0 - unhandledCropRate) =L=  maxLandChangeRate * ( suitableLandArea(location) -sum(crop_less_pasture, previousCropArea(crop_less_pasture, location))/*(1.0 - unhandledCropRate) - previousCropArea('pasture', location) );
-* PASTURE_DECREASE_CONSTRAINT(location) .. pastureDecrease(location) =L=  maxLandChangeRate * previousCropArea('pasture', location);
+ 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));
 
  COST_EQ .. total_cost =E=
          (
@@ -301,7 +290,7 @@ $gdxin
 
                SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
 
-               prematureDeforestationCost + SUM(location, totalConversionCost(location))
+               prematureDeforestationCost + SUM(location, totalConversionCost(location))  + A * 1000000
          );
 
  MODEL LAND_USE /ALL/ ;
@@ -309,7 +298,8 @@ $gdxin
  irrigI.L(crop, location) = previousIrrigIntensity(crop, location);
  otherIntensity.L(crop, location) = previousOtherIntensity(crop, location);
  cropArea.L(crop, location) = previousCropArea(crop, location);
- landCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, location);
+ landCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, location) - minimumLandCoverArea(land_cover, location);
+ restrictedLandCoverArea.L(land_cover, location) = minimumLandCoverArea(land_cover, location);
  ruminantFeed.L(feed_crop) = previousRuminantFeed(feed_crop);
  monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop);
  importAmount.L(all_types) = previousImportAmount(all_types);
@@ -329,6 +319,7 @@ $gdxin
  parameter netImportCost(all_types);
  parameter feedCostRate(feed_crop);
  parameter productionShock(all_types);
+ parameter reservedAreaDeforested(forest, location);
  scalar netCarbonFlux;
  scalar totalTimberProduction;
 
@@ -351,6 +342,7 @@ $gdxin
  netImportCost(import_crop) = importAmount.l(import_crop) * importPrices(import_crop) - exportAmount.l(import_crop) * exportPrices(import_crop);
  netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
 
+ reservedAreaDeforested(forest, location) = sum((land_cover)$[not sameAs(forest, land_cover)], restrictedLandCoverChange.L(forest, land_cover, location));
  netCarbonFlux = SUM(location, carbonFlux.L(location));
  totalTimberProduction = SUM(location, woodHarvest.L(location) + timberHarvest.L(location));
 

src/ac/ed/lurg/ModelConfig.java

@@ -249,7 +249,7 @@ public class ModelConfig {
 	// Wood/carbon data
 	public static final String FOREST_DIR = SPATIAL_DATA_DIR + File.separator + "forestry";
 	//public static final String CARBON_FLUX_FILE = FOREST_DIR + File.separator + "carbon_flux_";
-	//public static final String WOOD_YIELD_FILE = FOREST_DIR + File.separator + "wood_yield_";
+	public static final String WOOD_YIELD_FILENAME = "wood_yield.out";
 
 	// Output
 	public static final String LAND_COVER_OUTPUT_FILE = OUTPUT_DIR + File.separator + "lc.txt";

src/ac/ed/lurg/ModelMain.java

@@ -678,7 +678,7 @@ public class ModelMain {
 	private WoodYieldRasterSet getWoodYieldData(Timestep timestep) {
 		WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
 		LogWriter.println("Reading wood yield data");
-		new WoodYieldReader(woodYieldData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + "wood_yield.out");
+		new WoodYieldReader(woodYieldData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + ModelConfig.WOOD_YIELD_FILENAME);
 		return woodYieldData;
 	}
 

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

@@ -640,15 +640,15 @@ public class GamsLocationOptimiser {
 		
 			
 		// reserved Forest which was deforested
-		GAMSVariable varReservedDeforested = outDB.getVariable("reservedAreaDeforested");
+		GAMSParameter varReservedDeforested = outDB.getParameter("reservedAreaDeforested");
 		DoubleMap<Integer, LandCoverType, Double> reservedDeforested = new DoubleMap<Integer, LandCoverType, Double>();
 		
-		for (GAMSVariableRecord rec : varReservedDeforested) {
+		for (GAMSParameterRecord rec : varReservedDeforested) {
 			String landCover = rec.getKeys()[0];
 			String locationName = rec.getKeys()[1];
 			int locId = Integer.parseInt(locationName);
 			
-			double deforestedArea = rec.getLevel();
+			double deforestedArea = rec.getValue();
 			
 			if (deforestedArea > 0.00001) {
 				reservedDeforested.put(locId, LandCoverType.getForName(landCover), deforestedArea);

src/ac/ed/lurg/landuse/WoodYieldItem.java

@@ -8,6 +8,7 @@ import ac.sac.raster.RasterItem;
 
 public class WoodYieldItem implements RasterItem {
 	Map<LandCoverType, Map<LandCoverType, Double>> woodYields = new HashMap<LandCoverType, Map<LandCoverType, Double>>(); 
+	Map<LandCoverType, Map<LandCoverType, Double>> timberYields = new HashMap<LandCoverType, Map<LandCoverType, Double>>(); 
 	
 	public void setWoodYield(LandCoverType previousLandCover, LandCoverType newLandCover, double woodYield) {
 		if (woodYields.containsKey(previousLandCover)) {
@@ -19,10 +20,24 @@ public class WoodYieldItem implements RasterItem {
 		}
 	}
 	
+	public void setTimberYield(LandCoverType previousLandCover, LandCoverType newLandCover, double woodYield) {
+		if (timberYields.containsKey(previousLandCover)) {
+			timberYields.get(previousLandCover).put(newLandCover, woodYield);
+		} else {
+			Map<LandCoverType, Double> cfMap = new HashMap<LandCoverType, Double>();
+			cfMap.put(newLandCover, woodYield);
+			timberYields.put(previousLandCover, cfMap);
+		}
+	}
+	
 	public double getWoodYield(LandCoverType previousLandCover, LandCoverType newLandCover) {
 		return woodYields.get(previousLandCover).get(newLandCover);
 	}
 	
+	public double getTimberYield(LandCoverType previousLandCover, LandCoverType newLandCover) {
+		return timberYields.get(previousLandCover).get(newLandCover);
+	}
+	
 	public boolean checkForKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
 		if (woodYields.containsKey(previousLandCover)) {
 			if (woodYields.get(previousLandCover).containsKey(newLandCover)) {