Fixed ruminant overproduction bug.

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);

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

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);

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");

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);
 	}
 }

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);

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;