From e96ed8f233f819d78601061dc1661f68bbe65e6f Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Thu, 29 Apr 2021 20:30:07 +0100
Subject: [PATCH] Refactored GAMS code.
---
GAMS/IntExtOpt.gms | 96 +++++++++----------
src/ac/ed/lurg/ModelConfig.java | 2 +-
src/ac/ed/lurg/ModelMain.java | 2 +-
.../country/gams/GamsLocationOptimiser.java | 6 +-
src/ac/ed/lurg/landuse/WoodYieldItem.java | 15 +++
5 files changed, 64 insertions(+), 57 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index caecc0dd..883e6142 100644
--- a/GAMS/IntExtOpt.gms
+++ b/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));
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 7009cb27..59f788e7 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/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";
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index ca92fa0f..1aa7294b 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/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;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index ac22ad3f..078e720b 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/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);
diff --git a/src/ac/ed/lurg/landuse/WoodYieldItem.java b/src/ac/ed/lurg/landuse/WoodYieldItem.java
index 3da28759..87214369 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldItem.java
+++ b/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)) {
--
GitLab