From 8d8be0e19933669a44308995f0121e751ee1b47a Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Fri, 12 Mar 2021 18:16:20 +0000
Subject: [PATCH] Bug fixes. Code cleaned up.
---
GAMS/IntExtOpt.gms | 15 ++++++--
src/ac/ed/lurg/ModelMain.java | 5 ++-
src/ac/ed/lurg/country/CountryAgent.java | 4 ++
.../country/gams/GamsLocationOptimiser.java | 37 ++++++-------------
.../ed/lurg/country/gams/GamsRasterInput.java | 1 -
.../country/gams/GamsRasterOptimiser.java | 10 -----
src/ac/ed/lurg/landuse/LandUseItem.java | 20 +---------
src/ac/ed/lurg/landuse/WoodYieldReader.java | 1 -
src/ac/ed/lurg/types/LandCoverType.java | 26 ++++++++-----
src/ac/ed/lurg/utils/DoubleMap.java | 4 +-
src/ac/ed/lurg/utils/NestedMap.java | 3 --
src/ac/ed/lurg/utils/TripleMap.java | 3 +-
12 files changed, 49 insertions(+), 80 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 7cad291a..e4911ffe 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -120,7 +120,7 @@ $gdxin
*PARAMETER lcTransitionCost(lc_type_previous, lc_type, location) cost of land transitions including carbon cost and conversion cost;
*lcTransitionCost(lc_type_previous, lc_type, location) = carbonFluxRate(lc_type_previous, lc_type, location) * carbonPrice + agriExpansionCost(location);
- SCALAR forestExpansionMaxRate / 0.01 /;
+ SCALAR forestExpansionMaxRate / 0.05 /;
VARIABLES
cropArea(crop, location) total area for crops - Mha
@@ -144,12 +144,16 @@ $gdxin
woodHarvest(location) total wood harvested
carbonFlux(location) total carbon flux
totalFeedDM
+* A artificial variable "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,
agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease, totalFeedDM,
landCoverArea, landCoverChange, woodHarvest;
+* POSITIVE VARIABLE A;
+
EQUATIONS
UNIT_COST_EQ(crop, location) cost per area - $1000 per ha or $billion per Mha
YIELD_EQ(crop, location) yield given chosen intensity - tonnes per hectare
@@ -223,7 +227,10 @@ $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_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate) + cropArea('pasture', location);
+
+* Needs inequality due to floating point errors
+ TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(land_cover, landCoverArea(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);
@@ -254,9 +261,9 @@ $gdxin
FOREST_EXPANSION_CONSTRAINT(forest, location) .. landCoverArea(forest, location) - previousLandCoverArea(forest, location) =L= suitableLandArea(location) * forestExpansionMaxRate;
-* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverArea(land_cover, location) =G= minimumLandCoverArea(land_cover, location);
* minimum of minimumLandCoverArea and previousLandCoverArea corrects for floating point errors
- MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =G= min(minimumLandCoverArea(land_cover, location), previousLandCoverArea(land_cover, location));
+* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =G= min(minimumLandCoverArea(land_cover, location), previousLandCoverArea(land_cover, location));
+ MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =G= minimumLandCoverArea(land_cover, location);
WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYield(land_cover_before, land_cover_after, location));
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 66294c60..4bc0ff98 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -141,8 +141,7 @@ public class ModelMain {
double gen2Increase = (gen2EcDDemand>previousGen2EcDDemand) ? gen2EcDDemand - previousGen2EcDDemand : 0.0;
CarbonFluxRasterSet currentCarbonFluxData = getCarbonFluxData(timestep);
- LogWriter.println("Carbon flux data read");
- LogWriter.println("First line: " + currentCarbonFluxData.get(1, 1).getCarbonFluxes());
+
WoodYieldRasterSet currentWoodYieldData = getWoodYieldData(timestep);
countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, currentWoodYieldData);
@@ -666,6 +665,7 @@ public class ModelMain {
private CarbonFluxRasterSet getCarbonFluxData(Timestep timestep) {
CarbonFluxRasterSet carbonFluxData = new CarbonFluxRasterSet(desiredProjection);
new CarbonFluxReader(carbonFluxData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + "cflux_net.out");
+ LogWriter.println("Carbon flux data read");
return carbonFluxData;
}
@@ -673,6 +673,7 @@ public class ModelMain {
private WoodYieldRasterSet getWoodYieldData(Timestep timestep) {
WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
new WoodYieldReader(woodYieldData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + "wood_yield.out");
+ LogWriter.println("Wood yield data read");
return woodYieldData;
}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 0ab01089..6c55218a 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -160,6 +160,10 @@ public class CountryAgent extends AbstractCountryAgent {
FileSystems.getDefault().getPath(ModelConfig.TEMP_DIR + File.separator + "_gams_java_gdb2.gdx"),
FileSystems.getDefault().getPath(ModelConfig.TEMP_DIR + File.separator + country.getName().replaceAll("\\s+","") + currentTimestep.getYear() + ext + ".gdx"),
StandardCopyOption.REPLACE_EXISTING);
+ Files.copy(
+ FileSystems.getDefault().getPath(ModelConfig.TEMP_DIR + File.separator + "_gams_java_gjo1.lst"),
+ FileSystems.getDefault().getPath(ModelConfig.TEMP_DIR + File.separator + country.getName().replaceAll("\\s+","") + currentTimestep.getYear() + ext + "Listing" + ".lst"),
+ StandardCopyOption.REPLACE_EXISTING);
} catch (IOException e) {
LogWriter.print(e);
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 8fa74210..6844c309 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -31,7 +31,6 @@ import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.WoodYieldItem;
import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.GamsLandCoverType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.Parameter;
import ac.ed.lurg.types.YieldType;
@@ -126,7 +125,7 @@ public class GamsLocationOptimiser {
double suitableLand = landUseItem.getSuitableArea();
totalSuitable += suitableLand;
if (DEBUG) LogWriter.println(String.format(" %d %15s,\t %.3f", locationId, "suitableLand", suitableLand));
- setGamsParamValueTruncate(landP.addRecord(locString), suitableLand, 3);
+ setGamsParamValueTruncate(landP.addRecord(locString), suitableLand, 5);
double agriExpansionCost = ModelConfig.AGRI_EXPANSION_COST_BASE + landUseItem.getForestManagedFraction() * ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST;
if (DEBUG) LogWriter.println(String.format(" %d %15s,\t %.3f", locationId, "agriExpansionCost", agriExpansionCost));
@@ -176,33 +175,13 @@ public class GamsLocationOptimiser {
}
// Previous land covers
- // TODO Increase number of decimal places?
for (LandCoverType lc : LandCoverType.getConvertibleTypes()) {
Vector<String> v = new Vector<String>();
v.add(lc.getName());
v.add(locString);
setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(lc), 5);
}
-
- /*
- Vector<String> v = new Vector<String>();
-
- v.add("cropland");
- v.add(locString);
- setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(LandCoverType.CROPLAND), 3);
-
- v.set(0, "pasture");
- setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(LandCoverType.PASTURE), 3);
-
- v.set(0, "timberForest");
- setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST), 3);
-
- v.set(0, "carbonForest");
- setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(LandCoverType.CARBON_FOREST), 3);
-
- v.set(0, "natural");
- setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(LandCoverType.OTHER_NATURAL) + landUseItem.getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST), 3);
- */
+
}
if (DEBUG) LogWriter.println(String.format(" Total agricultural %.1f,\t suitable %.1f", totalAgriLand, totalSuitable));
@@ -417,9 +396,14 @@ public class GamsLocationOptimiser {
for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
Vector<String> v = new Vector<String>();
- v.add(coverMap.getKey().getName());
+ LandCoverType landCover = coverMap.getKey();
+ v.add(landCover.getName());
v.add(locString);
- setGamsParamValueTruncate(minimumLandCover.addRecord(v), coverMap.getValue(), 5);
+ double minCover = coverMap.getValue();
+ double previousCover = inputData.getPreviousLandUse().get(locationId).getUnprotectedLandCoverArea(landCover);
+ // Correction for floating point errors
+ double minCoverCorrected = Math.min(minCover, previousCover);
+ setGamsParamValueTruncate(minimumLandCover.addRecord(v), minCoverCorrected, 5);
}
}
@@ -649,13 +633,14 @@ public class GamsLocationOptimiser {
return results;
}
+ /*
private double getPrevUnmanagedForestProp(Integer locId) {
double prevOtherNaturalArea = inputData.getPreviousLandUse().get(locId).getUnprotectedLandCoverArea(LandCoverType.OTHER_NATURAL);
double prevUnmanagedForestArea = inputData.getPreviousLandUse().get(locId).getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST);
double prevUnmanagedForestProp = (prevOtherNaturalArea + prevUnmanagedForestArea) != 0 ? prevUnmanagedForestArea / (prevUnmanagedForestArea + prevOtherNaturalArea) : 0;
return prevUnmanagedForestProp;
}
-
+ */
private double getParmValue(GAMSParameter aParm, String itemName) {
try {
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index 1f655776..759107b3 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -5,7 +5,6 @@ import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.WoodYieldItem;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.DoubleMap;
-import ac.ed.lurg.utils.TripleMap;
import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.yield.YieldRaster;
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index d732bf50..74fe24d2 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -12,7 +12,6 @@ import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.WoodYieldItem;
import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.GamsLandCoverType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.YieldType;
import ac.ed.lurg.utils.DoubleMap;
@@ -219,11 +218,6 @@ public class GamsRasterOptimiser {
}
- if (locId == 31) {
- int foo = 1;
- }
-
-
// Allocate land cover change to grid cells
for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap: landCoverChange.getMap().entrySet()) {
LandCoverType fromLC = fromMap.getKey();
@@ -416,10 +410,6 @@ public class GamsRasterOptimiser {
CarbonFluxItem aggCFlux = aggregatedCarbonFluxes.lazyGet(clusterId);
WoodYieldItem aggWYield = aggregatedWoodYields.lazyGet(clusterId);
- if (clusterId == 1) {
- int foo = 1;
- }
-
landUseItem.updateSuitableArea(rasterInputData.getTimestep().getYear());
double suitableAreaThisTime = landUseItem.getSuitableArea();
double suitableAreaSoFar = aggLandUse.getSuitableArea();
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index b6e60f82..07541d8f 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -1,17 +1,13 @@
package ac.ed.lurg.landuse;
import java.io.Serializable;
-import java.util.ArrayList;
import java.util.Collection;
import java.util.HashMap;
import java.util.Map;
-import java.util.Map.Entry;
-
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.types.CropToDouble;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.utils.Interpolator;
import ac.sac.raster.InterpolatingRasterItem;
@@ -22,7 +18,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>(); //Mha
private Map<LandCoverType, Double> unprotectedAreas = new HashMap<LandCoverType, Double>(); //Mha
- //private DoubleMap<LandCoverType, Integer, Double> minimumLandCover = new DoubleMap<LandCoverType, Integer, Double>(); //Minimum land cover by year of expiry (Mha)
private double protectedArea; //protected area in Mha
private double unavailableArea; //area unavailable due to altitude etc
private double suitableArea; //Mha
@@ -291,7 +286,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
public double getTotalUnprotectedNatural() {
double unprotectedNatural = 0;
for (LandCoverType landType : LandCoverType.values()) {
- if (landType.isProtectable()) {
+ if (landType.isNatural()) {
unprotectedNatural += unprotectedAreas.get(landType);
}
}
@@ -436,19 +431,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return total;
}
- /*
- public double getMinimumLandCover(LandCoverType lcType, Integer year) {
- double area = 0;
- Map<Integer, Double> yearMap = minimumLandCover.getInnerMap(lcType);
- for (Entry<Integer, Double> entry : yearMap.entrySet()) {
- if (entry.getKey() > year) {
- area += entry.getValue();
- }
- }
- return area;
- }
- */
-
@Override
public String toString() {
return "LandUseItem: [landCoverAreas=" + landCoverAreas + ", protectedArea=" + protectedArea + ", unavailableArea=" + unavailableArea + "]";
diff --git a/src/ac/ed/lurg/landuse/WoodYieldReader.java b/src/ac/ed/lurg/landuse/WoodYieldReader.java
index eca499a9..76f92aaf 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldReader.java
+++ b/src/ac/ed/lurg/landuse/WoodYieldReader.java
@@ -2,7 +2,6 @@ package ac.ed.lurg.landuse;
import java.util.Map;
-import ac.ed.lurg.types.GamsLandCoverType;
import ac.ed.lurg.types.LandCoverType;
import ac.sac.raster.AbstractTabularRasterReader;
import ac.sac.raster.RasterKey;
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index 93d1fd0a..2784a4f6 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -9,24 +9,26 @@ import ac.ed.lurg.utils.LogWriter;
public enum LandCoverType {
- TIMBER_FOREST("timberForest", false, true, true),
- CARBON_FOREST("carbonForest", false, true, true),
- UNMANAGED_FOREST ("unmanagedForest", true, true, true),
- OTHER_NATURAL ("otherNatural", true, true, false),
- CROPLAND ("cropland", false, true, false),
- PASTURE ("pasture", false, true, false),
- BARREN ("barren", false, false, false),
- URBAN("urban", false, false, false);
+ TIMBER_FOREST("timberForest", false, true, true, true),
+ CARBON_FOREST("carbonForest", false, true, true, true),
+ UNMANAGED_FOREST ("unmanagedForest", true, true, true, true),
+ OTHER_NATURAL ("otherNatural", true, true, true, false),
+ CROPLAND ("cropland", false, true, false, false),
+ PASTURE ("pasture", false, true, false, false),
+ BARREN ("barren", false, false, false, false),
+ URBAN("urban", false, false, false, false);
private String name;
private boolean isProtectable;
private boolean isConvertible;
+ private boolean isNatural;
private boolean isForest;
- LandCoverType(String name, boolean isProtectable, boolean isConvertible, boolean isForest) {
+ LandCoverType(String name, boolean isProtectable, boolean isConvertible, boolean isNatural, boolean isForest) {
this.name = name;
this.isProtectable = isProtectable;
this.isConvertible = isConvertible;
+ this.isNatural = isNatural;
this.isForest = isForest;
}
@@ -38,6 +40,10 @@ public enum LandCoverType {
return isProtectable;
}
+ public boolean isNatural() {
+ return isNatural;
+ }
+
private static final Map<String, LandCoverType> nameCache = new HashMap<String, LandCoverType>();
static {
for (LandCoverType c : values()) {
@@ -80,6 +86,7 @@ public enum LandCoverType {
}
+ /*
public static Collection<LandCoverType> getConvertibleNonForestTypes() {
Collection<LandCoverType> convNonforestTypes = new HashSet<LandCoverType>();
@@ -91,5 +98,6 @@ public enum LandCoverType {
return convNonforestTypes;
}
+ */
}
diff --git a/src/ac/ed/lurg/utils/DoubleMap.java b/src/ac/ed/lurg/utils/DoubleMap.java
index 14f7fdcf..8a50e7c0 100644
--- a/src/ac/ed/lurg/utils/DoubleMap.java
+++ b/src/ac/ed/lurg/utils/DoubleMap.java
@@ -6,7 +6,7 @@ import java.util.Map;
public class DoubleMap<K, L, V> {
// Stores a double mapped by two consecutive keys
- Map<K, Map<L, Double>> theMap;
+ private Map<K, Map<L, Double>> theMap;
public DoubleMap() {
this.theMap = new HashMap<K, Map<L, Double>>();
@@ -50,7 +50,5 @@ public class DoubleMap<K, L, V> {
} else {
put(key1, key2, value);
}
-
-
}
}
diff --git a/src/ac/ed/lurg/utils/NestedMap.java b/src/ac/ed/lurg/utils/NestedMap.java
index 397f686f..4fc2bebc 100644
--- a/src/ac/ed/lurg/utils/NestedMap.java
+++ b/src/ac/ed/lurg/utils/NestedMap.java
@@ -1,9 +1,6 @@
package ac.ed.lurg.utils;
-import java.util.Collections;
import java.util.HashMap;
-import java.util.Map;
-import java.util.List;
public class NestedMap<K, V> {
diff --git a/src/ac/ed/lurg/utils/TripleMap.java b/src/ac/ed/lurg/utils/TripleMap.java
index 11c6cf37..aa9e32dd 100644
--- a/src/ac/ed/lurg/utils/TripleMap.java
+++ b/src/ac/ed/lurg/utils/TripleMap.java
@@ -5,7 +5,7 @@ import java.util.Map;
public class TripleMap<K, L, M, V> {
- Map<K, DoubleMap<L, M, V >> theMap = new HashMap<K, DoubleMap<L, M, V>>();
+ private Map<K, DoubleMap<L, M, V >> theMap = new HashMap<K, DoubleMap<L, M, V>>();
public TripleMap() {
@@ -21,7 +21,6 @@ public class TripleMap<K, L, M, V> {
}
}
-
public Double get(K key1, L key2, M key3) {
if (theMap.containsKey(key1)) {
return theMap.get(key1).get(key2, key3);
--
GitLab