From c997e73fcbba19ef11446474f9d9687d74ee5683 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Wed, 16 Sep 2015 20:19:30 +0100
Subject: [PATCH] Corrections to GAMS energy calc

---
 GAMS/IntExtOpt.gms                            | 18 +++++----
 src/ac/ed/lurg/ModelConfig.java               |  2 +-
 .../lurg/country/CompositeCountryManager.java |  2 +-
 .../country/gams/GamsLocationOptimiser.java   | 10 ++++-
 src/ac/ed/lurg/demand/BaseConsumpManager.java |  2 +-
 src/ac/ed/lurg/demand/DemandManager.java      | 38 +++++++++----------
 6 files changed, 41 insertions(+), 31 deletions(-)

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 750f4a0e..db113c5c 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -62,7 +62,7 @@ $gdxin
               oilcrops          0.7 
               pulses            0.5
               starchyRoots      4.0
-              pasture           0.1  / ; 
+              pasture           0.0  / ; 
                
  VARIABLES
        area(crop, location)               total area for each crop - Mha
@@ -163,12 +163,16 @@ $gdxin
  PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= area('pasture', location) - previousArea('pasture', location);
  PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(area('pasture', location) - previousArea('pasture', location));
   
- ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location))     
-            + (sum(location, agriLandExpansion(location)) * landChangeEnergy)
-            + (sum(location, (0.2 * SQR((cropIncrease(location) + cropDecrease(location)) / suitableLandArea(location)) + (cropIncrease(location) + cropDecrease(location)) * 0.1) * (cropIncrease(location) + cropDecrease(location)) * landChangeEnergy))
-            + (sum(location, (0.2 * SQR((pastureIncrease(location) + pastureDecrease(location)) / suitableLandArea(location)) + (pastureIncrease(location) + pastureDecrease(location)) * 0.1) * (pastureIncrease(location) + pastureDecrease(location)) * landChangeEnergy))
-            + sum(import_crop, (importAmount(import_crop)) * worldImportPrices(import_crop))
-            - sum(import_crop, (exportAmount(import_crop)) * worldExportPrices(import_crop))) / 1000000;          
+ ENERGY_EQ .. energy =E= 
+         ( 
+              SUM((crop, location), area(crop, location) * unitEnergy(crop, location)) +
+              sum(location, agriLandExpansion(location) * landChangeEnergy) +
+              sum(location, (SQR(cropIncrease(location) + cropDecrease(location)) / (suitableLandArea(location) * 0.01 + sum(crop_less_pasture, previousArea(crop_less_pasture, location)))
+                     + (cropIncrease(location) + cropDecrease(location)) * 0.1) * landChangeEnergy) + 
+              sum(location, (SQR(pastureIncrease(location) + pastureDecrease(location)) / (suitableLandArea(location) * 0.01 + previousArea('pasture', location))
+                     + (pastureIncrease(location) + pastureDecrease(location)) * 0.1) * landChangeEnergy) +
+              sum(import_crop, importAmount(import_crop) * worldImportPrices(import_crop) - exportAmount(import_crop) * worldExportPrices(import_crop)) 
+         ) / 1000000;          
  
  MODEL LAND_USE /ALL/ ;
  
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 8c6085b0..c507777c 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -149,7 +149,7 @@ public class ModelConfig {
 	public static final double POPULATION_AGGREG_LIMIT = getDoubleProperty("POPULATION_AGGREG_LIMIT", 40.0);  // in millions, smaller countries are aggregated on a regional basis
 	
 	public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 2.0);
-	public static final double FERTILISER_MAX_COST = getDoubleProperty("FERTILISER_MAX_COST", 1.0);
+	public static final double FERTILISER_MAX_COST = getDoubleProperty("FERTILISER_MAX_COST", 2.5);
 	public static final double TRANSPORT_LOSSES = getDoubleProperty("TRANSPORT_LOSSES", 0.15);  // in international trade
 	public static final double TRADE_BARRIER_FACTOR = getDoubleProperty("TRADE_BARRIER_FACTOR", 1.2);  // price factor in international trade, transport cost and real trade barriers
 
diff --git a/src/ac/ed/lurg/country/CompositeCountryManager.java b/src/ac/ed/lurg/country/CompositeCountryManager.java
index 34a46119..64352124 100644
--- a/src/ac/ed/lurg/country/CompositeCountryManager.java
+++ b/src/ac/ed/lurg/country/CompositeCountryManager.java
@@ -25,7 +25,7 @@ public class CompositeCountryManager {
 	private void populate(BaseConsumpManager baseConsumpManager) {
 		mapFromSingleCountry = new HashMap<SingleCountry, CompositeCountry>();
 		
-		HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("Bangladesh", "Portugal", "Haiti"));
+		HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("Bangladesh", "Portugal", "Haiti", "Democratic Republic of the Congo"));
 		
 		for (SingleCountry c : baseConsumpManager.getAllCountries()) {
 			CompositeCountry cc;
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index c03b8ade..06ee659b 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -170,10 +170,16 @@ public class GamsLocationOptimiser {
 			LogWriter.println("\nImport-export, min net imports, max net imports,  import price,   export price");
 			for (CropType crop : CropType.getImportedTypes()) {
 				ImportExportConstraint iec = countryInput.getImportConstraints().get(crop);
-				GlobalPrice gp = countryInput.getWorldPrices().get(crop);
+				double minNetImport = 0, maxNetImport = 0;
+				if (iec != null) {
+					minNetImport = iec.getMinNetImport();
+					maxNetImport = iec.getMaxNetImport();
+				}
 				
+				GlobalPrice gp = countryInput.getWorldPrices().get(crop);
+								
 				LogWriter.println(String.format("     %15s, \t %5.1f, \t %5.1f, \t %5.3f, \t %5.3f", 
-						crop.getGamsName(), iec.getMinNetImport(), iec.getMaxNetImport(), gp.getImportPrice(), gp.getExportPrice()));				
+						crop.getGamsName(), minNetImport, maxNetImport, gp.getImportPrice(), gp.getExportPrice()));
 			}
 		}
 		
diff --git a/src/ac/ed/lurg/demand/BaseConsumpManager.java b/src/ac/ed/lurg/demand/BaseConsumpManager.java
index 8af8c4c7..e79d7123 100644
--- a/src/ac/ed/lurg/demand/BaseConsumpManager.java
+++ b/src/ac/ed/lurg/demand/BaseConsumpManager.java
@@ -71,7 +71,7 @@ public class BaseConsumpManager {
 
 		if (commodityMap == null || !commodityMap.containsKey(commodity)) {
 			LogWriter.printlnError("BaseConsumpManager: can't get value for " + commodity + ", " + country);
-			return Double.NaN;
+			return 0.0;
 		}
 		
 		Double d = commodityMap.get(commodity);
diff --git a/src/ac/ed/lurg/demand/DemandManager.java b/src/ac/ed/lurg/demand/DemandManager.java
index e666003d..4e6c8316 100644
--- a/src/ac/ed/lurg/demand/DemandManager.java
+++ b/src/ac/ed/lurg/demand/DemandManager.java
@@ -18,7 +18,7 @@ public class DemandManager {
 	private BaseConsumpManager baseConsumpManager;
 	private CompositeCountryManager compositeCountryManager;
 	private Map<SingleCountry, Map<CommodityType, Double>> bioenergyBaseDemand;
-	
+
 	private ModelFitType fitType;
 	private String ssp_scenario;
 
@@ -31,46 +31,46 @@ public class DemandManager {
 		sspManager = new SspManager();
 		bioenergyBaseDemand = new BioenergyDemand().getBioenergyData();
 	}
-	
+
 	public Map<CommodityType, Double> getDemand(CompositeCountry cc, int year) {
-		
+
 		if (ModelConfig.KEEP_DEMAND_FIXED)
 			year = ModelConfig.BASE_YEAR;
-		
+
 		Map<CommodityType, Double> demandMap = new HashMap<CommodityType, Double>();
 
 		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
-		SspData baseSspData = sspManager.get(ssp_scenario, ModelConfig.BASE_YEAR, c);
-		
+			SspData baseSspData = sspManager.get(ssp_scenario, ModelConfig.BASE_YEAR, c);
+
 			SspData sd = sspManager.get(ssp_scenario, year, c);
-			
+
 			LogWriter.println("Got ssp data for " + c.getCountryName() + " of " + sd);
-			
+
 			if (sd == null) {
 				LogWriter.printlnError("No ssp data for " + ssp_scenario + ", " + year + ", " + c.getCountryName());
 				continue;
 			}
-			
+
 			for (DemandCurve dc : demandCurveManager.get(fitType)) {
 				double cpc = dc.getConsumptionPc(baseSspData.getGdpPc(), baseConsumpManager.get(c, dc.getCommodityType()), sd.getGdpPc());
 				double d = cpc * sd.getPopulation();
-				
+
 				double bioenergy = getBioenergyDemand(c, year, dc.getCommodityType());
-		//		LogWriter.println(String.format("Country %s comm %s: %f",  country, dc.getCommodityType(), bioenergy));
-				
+				//		LogWriter.println(String.format("Country %s comm %s: %f",  country, dc.getCommodityType(), bioenergy));
+
 				double newDemand = (d + bioenergy) * (1.0 + ModelConfig.SEED_AND_WASTE_FRACTION);
 				Double prevDemand = demandMap.get(dc.getCommodityType());
-				
+
 				if (prevDemand != null)
 					newDemand += prevDemand.doubleValue();  // aggregate if composite country
-				
+
 				demandMap.put(dc.getCommodityType(), newDemand);
 			}
 		}
-		
+
 		return demandMap;
 	}
-	
+
 	private double getBioenergyDemand(SingleCountry country, int year, CommodityType commodity) {
 		// could adjust for year somehow, but not doing this yet
 		if (bioenergyBaseDemand != null && bioenergyBaseDemand.containsKey(country)) {
@@ -79,8 +79,8 @@ public class DemandManager {
 		}
 		return 0.0;
 	}
-	
-/*	public double getPopulation (SingleCountry country, int year) {
+
+	/*	public double getPopulation (SingleCountry country, int year) {
 		SspData sd = sspManager.get(ssp_scenario, year, country);
 		if (sd == null) {
 			LogWriter.printlnError("No ssp data for " + ssp_scenario + ", " + year + ", " + country.getCountryName());
@@ -93,5 +93,5 @@ public class DemandManager {
 	public String getModelAndScenarioDescription() {
 		return fitType + " " + ssp_scenario;
 	}
-	*/
+	 */
 }
-- 
GitLab