From c137553c1cc21443609b73e0d962f997a3754fd8 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Sat, 2 Apr 2022 21:05:23 +0100
Subject: [PATCH] Fixed ruminant overproduction bug.
---
GAMS/IntExtOpt.gms | 38 +++++++++++--------
src/ac/ed/lurg/ModelConfig.java | 6 +--
src/ac/ed/lurg/ModelMain.java | 4 +-
.../country/gams/GamsDemandOptimiser.java | 1 +
src/ac/ed/lurg/landuse/LandTileReader.java | 4 +-
src/ac/ed/lurg/landuse/WoodUsageReader.java | 4 +-
src/ac/ed/lurg/types/WoodType.java | 2 +-
7 files changed, 34 insertions(+), 25 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index eb9064de..306186cb 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -122,18 +122,19 @@ $gdxin
* minNetImport(import_crop) = min(minNetImport(import_crop), demand(import_crop));
VARIABLES
- cropArea(crop, location) total area for crops - Mha
+ unitCost(crop, location) cost per area for each crop - cost
fertI(crop, location) fertilizer intensity for each crop - factor between 0 and 1
irrigI(crop, location) irrigation intensity for each crop - factor between 0 and 1
otherIntensity(crop, location)
+ yield(crop, location) yield per area for each crop - t per ha
+ cropArea(crop, location) total area for crops - Mha
+ netSupply(all_types) crop supply
ruminantFeed(crop) amount of feed for ruminant animals - Mt
monogastricFeed(crop) amount of feed for monogatric animals - Mt
importAmount(all_types) imports of crops and meat - Mt
exportAmount(all_types) exports of crops and meat - Mt
- yield(crop, location) yield per area for each crop - t per ha
- unitCost(crop, location) cost per area for each crop - cost
+
totalFeedDM total feed dry matter
- supply(all_types) crop supply
landCoverArea(land_cover, location) land cover area in Mha
landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha
totalConversionCost(location) land cover conversion cost - $1000 per ha
@@ -163,12 +164,14 @@ $gdxin
MAX_OTHER_INTENSITY_CONSTRAINT(crop, location)
IRRIGATION_CONSTRAINT(location) constraint on water usage
- SUPPLY_CALC(crop)
- RUMINANTS_SUPPLY_CALC
- MONOGASTRICS_SUPPLY_CALC
+ NET_SUPPLY_EQ(crop) calc net supply for crops
+ CROP_DEMAND_CONSTRAINT(crop) satisfy demand for individual crops
+ TOTAL_CEREAL_DEMAND_CONSTRAINT satisfy demand for combined cereals
+ TOTAL_OIL_PULSE_DEMAND_CONSTRAINT satisfy demand for combined oilcrops and pulses
+ RUMINANT_DEMAND_CONSTRAINT satisfy demand for ruminant products
+ MONOGASTRICS_DEMAND_CONSTRAINT satisfy demand for monogastric products
TOTAL_NON_PASTURE_FEED_DM_CALC calc total feed dry matter not including pasture
FEED_MIX_CONSTRAINT(feed_crop_less_pasture) limit amount of feed for each feed crop
- SUPPLY_CONSTRAINT(import_crop)
MAX_NET_IMPORT_CONSTRAINT(import_crop) constraint on max net imports
MIN_NET_IMPORT_CONSTRAINT(import_crop) constraint on min net imports
@@ -221,16 +224,21 @@ $gdxin
*************** Commodity supply *************************
- SUPPLY_CALC(crop) .. supply(crop) =E= sum(location, cropArea(crop, location) * yield(crop, location)) * (1 - seedAndWasteRate(crop)) - ruminantFeed(crop) - monogastricFeed(crop);
+ NET_SUPPLY_EQ(crop) .. netSupply(crop) =E= sum(location, cropArea(crop, location) * yield(crop, location)) * (1 - seedAndWasteRate(crop)) - ruminantFeed(crop) - monogastricFeed(crop) + importAmount(crop) - exportAmount(crop);
+
+ CROP_DEMAND_CONSTRAINT(crop) .. netSupply(crop) =G= demand(crop);
+
+ TOTAL_CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, netSupply(cereal_crop)) =G= demand('cereals');
- RUMINANTS_SUPPLY_CALC .. supply('ruminants') =E= meatEfficency * sum(feed_crop, ruminantFeed(feed_crop) * cropDM(feed_crop)) * (1 - seedAndWasteRate('ruminants'));
+ TOTAL_OIL_PULSE_DEMAND_CONSTRAINT .. sum(oilpulse_crop, netSupply(oilpulse_crop)) =G= demand('oilcropspulses');
+
+ RUMINANT_DEMAND_CONSTRAINT .. meatEfficency * sum(feed_crop, ruminantFeed(feed_crop) * cropDM(feed_crop)) * (1 - seedAndWasteRate('ruminants')) =G= (demand('ruminants') - importAmount('ruminants') + exportAmount('ruminants'));
+
+ MONOGASTRICS_DEMAND_CONSTRAINT .. meatEfficency * sum(feed_crop_less_pasture, monogastricFeed(feed_crop_less_pasture) * cropDM(feed_crop_less_pasture)) * (1 - seedAndWasteRate('monogastrics')) =G= (demand('monogastrics') - importAmount('monogastrics') + exportAmount('monogastrics'));
- MONOGASTRICS_SUPPLY_CALC .. supply('monogastrics') =E= meatEfficency * sum(feed_crop_less_pasture, monogastricFeed(feed_crop_less_pasture) * cropDM(feed_crop_less_pasture)) * (1 - seedAndWasteRate('monogastrics'));
-
TOTAL_NON_PASTURE_FEED_DM_CALC .. totalFeedDM =E= sum(feed_crop_less_pasture, (ruminantFeed(feed_crop_less_pasture) + monogastricFeed(feed_crop_less_pasture)) * cropDM(feed_crop_less_pasture));
FEED_MIX_CONSTRAINT(feed_crop_less_pasture) .. (ruminantFeed(feed_crop_less_pasture) + monogastricFeed(feed_crop_less_pasture)) * cropDM(feed_crop_less_pasture) =L= totalFeedDM * 0.7;
-
- SUPPLY_CONSTRAINT(import_crop) .. supply(import_crop) =E= demand(import_crop) + exportAmount(import_crop) - importAmount(import_crop);
+
*************** Exports *******************************
@@ -313,7 +321,7 @@ $gdxin
SOLVE LAND_USE USING NLP MINIMIZING total_cost;
display previousCropArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, cropArea.L;
- display supply.l, demand, ruminantFeed.l, monogastricFeed.l, exportAmount.l;
+ display netSupply.l, demand, ruminantFeed.l, monogastricFeed.l, exportAmount.l;
* Calculate summary information used in Java process
parameter totalProd(all_types);
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 8c7d7116..8e274af0 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -462,9 +462,8 @@ public class ModelConfig {
public static final String TIMBER_DEMAND_FILE = getProperty("TIMBER_DEMAND_FILE", DATA_DIR + File.separator + TIMBER_DEMAND_FILENAME);
public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 1000.0); // tC-eq
- public static final double ROUNDWOOD_DEMAND_ELASTICITY = getDoubleProperty("ROUNDWOOD_DEMAND_ELASTICITY", -0.21956);
- public static final double FUELWOOD_DEMAND_ELASTICITY = getDoubleProperty("FUELDWOOD_DEMAND_ELASTICITY", -0.22787);
- public static final double FOREST_ESTABLISHMENT_COST = getDoubleProperty("FOREST_ESTABLISHMENT_COST", 2.0); // $1000/ha
+ public static final double IND_ROUNDWOOD_DEMAND_ELASTICITY = getDoubleProperty("IND_ROUNDWOOD_DEMAND_ELASTICITY", 0.3123881);
+ public static final double FUELWOOD_DEMAND_ELASTICITY = getDoubleProperty("FUELDWOOD_DEMAND_ELASTICITY", -0.3598551);
public static final String WOOD_NET_IMPORTS_FILENAME = getProperty("WOOD_NET_IMPORTS_FILENAME", "wood_net_imports.csv");
public static final String WOOD_NET_IMPORTS_FILE = getProperty("WOOD_NET_IMPORTS_FILE", DATA_DIR + File.separator + WOOD_NET_IMPORTS_FILENAME);
public static final double WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 3e-7); // m3 to MtC-eq p.16 [https://doi.org/10.5194/gmd-13-5425-2020]
@@ -473,6 +472,7 @@ public class ModelConfig {
public static final boolean CONVERSION_COSTS_FROM_FILE = getBooleanProperty("CONVERSION_COSTS_FROM_FILE", false);
public static final boolean IS_FORESTRY_ON = getBooleanProperty("IS_FORESTRY_ON", true);
public static final boolean IS_CARBON_ON = getBooleanProperty("IS_CARBON_ON", true);
+ public static final double FOREST_ESTABLISHMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_ESTABLISHMENT_COST", 2.0) : 0.0; // $1000/ha
public static final int CARBON_HORIZON = getIntProperty("CARBON_HORIZON", 30);
public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.05); //$1000/tC
public static final double INIT_CARBON_PRICE = IS_CARBON_ON ? getDoubleProperty("INIT_CARBON_PRICE", 0.02) : 0.0; // $1000/tC-eq
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 83dc5606..3a2f4a15 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -140,7 +140,7 @@ public class ModelMain {
private void doTimestep(Timestep timestep) {
System.gc();
LogWriter.println("Timestep: " + timestep.toString());
- LogWriter.println("Memory usage 1: " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()));
+ LogWriter.println("Memory usage 1: " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()) / (1024.0*1024.0*1024.0));
YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
getUpdateIrrigationData(timestep, yieldSurfaces); // updating currentIrrigationData
@@ -183,7 +183,7 @@ public class ModelMain {
}
internationalMarket.applyPriceShocks(timestep);
- LogWriter.println("Memory usage: 2" + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()));
+ LogWriter.println("Memory usage 2: " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()) / (1024.0*1024.0*1024.0));
// output results
outputTimestepResults(timestep);
diff --git a/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java b/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
index 798fa03d..febfca4c 100755
--- a/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
@@ -58,6 +58,7 @@ public class GamsDemandOptimiser {
opt.defines("gdx_parameters", new File(ModelConfig.DEMAND_PARAM_FILE).getAbsolutePath());
long startTime = System.currentTimeMillis();
+ LogWriter.println("Memory usage: GAMS Demand " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()) / (1024.0*1024.0*1024.0));
gamsJob.run(opt, dbs.toArray(new GAMSDatabase[dbs.size()]));
LogWriter.println("Took " + (System.currentTimeMillis() - startTime) + " ms to run");
diff --git a/src/ac/ed/lurg/landuse/LandTileReader.java b/src/ac/ed/lurg/landuse/LandTileReader.java
index 4a3b3aae..f1772633 100644
--- a/src/ac/ed/lurg/landuse/LandTileReader.java
+++ b/src/ac/ed/lurg/landuse/LandTileReader.java
@@ -8,7 +8,7 @@ import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterSet;
public class LandTileReader extends AbstractTabularRasterReader<LandCoverItem> {
- private static final int MIN_COLS = 6;
+ private static final int MIN_COLS = 7;
public LandTileReader(RasterSet<LandCoverItem> data) {
super("[ |\t]+", MIN_COLS, data);
@@ -21,6 +21,6 @@ public class LandTileReader extends AbstractTabularRasterReader<LandCoverItem> {
item.setLandUseAgeDist(LandCoverType.PASTURE, getValueForCol(rowValues, "PASTURE"), ageGroup);
item.setLandUseAgeDist(LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "FOREST"), ageGroup);
item.setLandUseAgeDist(LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "FOREST"), ageGroup);
- item.setLandUseAgeDist(LandCoverType.NATURAL, getValueForCol(rowValues, "FOREST"), ageGroup);
+ item.setLandUseAgeDist(LandCoverType.NATURAL, getValueForCol(rowValues, "NATURAL"), ageGroup);
}
}
diff --git a/src/ac/ed/lurg/landuse/WoodUsageReader.java b/src/ac/ed/lurg/landuse/WoodUsageReader.java
index 0a3e50d7..add6e103 100644
--- a/src/ac/ed/lurg/landuse/WoodUsageReader.java
+++ b/src/ac/ed/lurg/landuse/WoodUsageReader.java
@@ -53,8 +53,8 @@ public class WoodUsageReader {
countryName = tokens[COUNTRY_COL];
woodTypeName = tokens[WOOD_TYPE_COL];
- harvest = Double.parseDouble(tokens[HARVEST_COL]) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; // m3 to MtC-eq
- netImport = Double.parseDouble(tokens[NET_IMPORT_COL]) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; // m3 to MtC-eq
+ harvest = ModelConfig.IS_FORESTRY_ON ? Double.parseDouble(tokens[HARVEST_COL]) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR : 0; // m3 to MtC-eq
+ netImport = ModelConfig.IS_FORESTRY_ON ? Double.parseDouble(tokens[NET_IMPORT_COL]) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR : 0; // m3 to MtC-eq
SingleCountry country = CountryManager.getForName(countryName);
WoodType woodType = WoodType.getForName(woodTypeName);
diff --git a/src/ac/ed/lurg/types/WoodType.java b/src/ac/ed/lurg/types/WoodType.java
index de4af65f..a8407076 100644
--- a/src/ac/ed/lurg/types/WoodType.java
+++ b/src/ac/ed/lurg/types/WoodType.java
@@ -5,7 +5,7 @@ import java.util.Map;
import ac.ed.lurg.ModelConfig;
public enum WoodType {
- ROUNDWOOD("roundwood", ModelConfig.ROUNDWOOD_DEMAND_ELASTICITY),
+ ROUNDWOOD("roundwood", ModelConfig.IND_ROUNDWOOD_DEMAND_ELASTICITY),
FUELWOOD("fuelwood", ModelConfig.FUELWOOD_DEMAND_ELASTICITY);
private String name;
--
GitLab