From baffd107b2db3bda2453f736d39a3e5f116101d1 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Fri, 14 Aug 2015 17:54:35 +0100
Subject: [PATCH] Update global prices based on balance of imports and exports

---
 GAMS/IntExtOpt.gms                            |  8 +--
 src/ac/ed/lurg/ModelConfig.java               |  7 +-
 src/ac/ed/lurg/ModelMain.java                 | 68 +++++++++++--------
 src/ac/ed/lurg/country/CountryAgent.java      |  8 +--
 .../lurg/country/gams/GamsCountryInput.java   | 14 ++--
 .../country/gams/GamsLocationOptimiser.java   |  4 +-
 6 files changed, 61 insertions(+), 48 deletions(-)

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 2a0baeee..6a08c228 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -19,7 +19,7 @@
  PARAMETER fertParam(crop, location)         yield response to fertilizer parameter;
  PARAMETER irrigParam(crop, location)        yield response to irrigation parameter;
  PARAMETER demand(all_types)                 in t;
- PARAMETER worldInputEnergy(import_crop)     average input energy from world exports used to determine if we should import or export;
+ PARAMETER worldInputPrices(import_crop)     average input energy from world exports used to determine if we should import or export;
  PARAMETER maxNetImport(import_crop)         maximum net import for each crop based on world market;
  PARAMETER minNetImport(import_crop)         minimum net import for each crop based on world market;
  PARAMETER irrigCost(location)               irrigation cost;
@@ -35,7 +35,7 @@
 $gdxin %gdxincname% 
 $load location, suitableLandArea, previousArea, demand, landChangeEnergy
 $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
-$load fertParam, irrigParam, worldInputEnergy, maxNetImport, minNetImport
+$load fertParam, irrigParam, worldInputPrices, maxNetImport, minNetImport
 $load meatEfficency, minFeedRate, irrigCost
 $load cropAdj, fitToPreviousAreas
 $gdxin
@@ -157,7 +157,7 @@ $gdxin
             + (sum(location, agriLandExpansion(location)) * landChangeEnergy)
             + (sum(location, (cropIncrease(location) + cropDecrease(location)) * landChangeEnergy/8))
             + (sum(location, (pastureIncrease(location) + pastureDecrease(location)) * landChangeEnergy/8 ))
-            + sum(import_crop, (netImportAmount(import_crop)) * worldInputEnergy(import_crop))) / 1000000;          
+            + sum(import_crop, (netImportAmount(import_crop)) * worldInputPrices(import_crop))) / 1000000;          
  
  MODEL LAND_USE /ALL/ ;
  
@@ -217,7 +217,7 @@ $gdxin
  parameter feedEnergy(all_types);
  
  totalProdCost(crop) = sum(location, unitEnergy.l(crop, location) * area.l(crop, location));
- totalImportEnergy(import_crop) = netImportAmount.l(import_crop) * worldInputEnergy(import_crop);
+ totalImportEnergy(import_crop) = netImportAmount.l(import_crop) * worldInputPrices(import_crop);
  
  feedEnergy(feed_crop)$(netImportAmount.l(feed_crop) gt 0) = (totalProdCost(feed_crop) + totalImportEnergy(feed_crop)) * feedAmount.l(feed_crop) / (totalProd(feed_crop) + netImportAmount.l(feed_crop));
  feedEnergy(feed_crop)$(netImportAmount.l(feed_crop) le 0 and totalProd(feed_crop) gt 0) = totalProdCost(feed_crop) * feedAmount.l(feed_crop) / totalProd(feed_crop);
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index db97a623..594d9e2d 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -118,8 +118,8 @@ public class ModelConfig {
 
 	
 	public static final int START_TIMESTEP = getIntProperty("START_TIMESTEP", 0);
-	public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 1);
-	public static final int TIMESTEP_SIZE = getIntProperty("END_TIMESTEP", 40);
+	public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 4);
+	public static final int TIMESTEP_SIZE = getIntProperty("TIMESTEP_SIZE", 10);
 	public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2010);
 	
 	public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", 0.1);
@@ -138,4 +138,7 @@ public class ModelConfig {
 	public static final double LAND_CHANGE_COST = getDoubleProperty("LAND_CHANGE_COST", 2.0);
 	public static final double MIN_FEED_RATE = getDoubleProperty("MIN_FEED_RATE", 0.15);
 	public static final double SEED_AND_WASTE_FRACTION = getDoubleProperty("SEED_AND_WASTE_FRACTION", 0.1);
+	
+	public static final double MARKET_LAMBA = getDoubleProperty("MARKET_LAMBA", 0.3); // controls market price adjustment rate
+
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index a6ff0124..546584a2 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -45,7 +45,7 @@ public class ModelMain {
 	private Collection<CountryAgent> countryAgents;
 	private Map<Country, List<RasterKey>> countryToKeysMap;
 	private DemandManager demandManager;
-	private Map<CropType, Double> prevWorldInputCost;
+	private Map<CropType, Double> prevWorldPrices;
 	private RasterHeaderDetails desiredProjection;
 
 	public static void main(String[] args)  {
@@ -62,9 +62,9 @@ public class ModelMain {
 		countryAgents = createCountryAgents();
 
 		// in first timestep we don't have this info, but ok as constrained to import/export specified amount
-		prevWorldInputCost = new HashMap<CropType, Double>();  
+		prevWorldPrices = new HashMap<CropType, Double>();  
 		for (CropType c : CropType.getImportedTypes()) 
-			prevWorldInputCost.put(c, 1.0);
+			prevWorldPrices.put(c, 1.0);
 	}
 
 	/* run the model */
@@ -90,11 +90,8 @@ public class ModelMain {
 		//		YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-50.0, -4.0);
 		//		LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.CEREALS)  + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
 
-
-		CropToDoubleMap totalQuantity = new CropToDoubleMap();
-		CropToDoubleMap totalWorldInputCost = new CropToDoubleMap();
-		//	CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
-		//	CropToDoubleMap totalExportCommodities = new CropToDoubleMap();
+		CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
+		CropToDoubleMap totalExportCommodities = new CropToDoubleMap();
 
 		RasterSet<IntensitiesItem> globalIntensityRaster = new RasterSet<IntensitiesItem>(desiredProjection);
 		RasterSet<AreasItem> globalCropAreaRaster = new RasterSet<AreasItem>(desiredProjection);
@@ -107,7 +104,7 @@ public class ModelMain {
 			YieldRaster countryYieldSurfaces = yieldSurfaces.getSubsetRasterForKeys(countryKeys);
 						
 			// do the optimization
-			GamsRasterOutput result = ca.determineProduction(timestep, countryYieldSurfaces, prevWorldInputCost);
+			GamsRasterOutput result = ca.determineProduction(timestep, countryYieldSurfaces, prevWorldPrices);
 
 			if (result == null) {
 				LogWriter.printlnError("No results for " + ca.getCountry());
@@ -124,34 +121,47 @@ public class ModelMain {
 
 			for (Entry<CropType, CropUsageData> entry : cropUsage.entrySet()) {
 				CropType c = entry.getKey();
-				CropUsageData d = entry.getValue();
+				double countryNetImports = entry.getValue().getNetImports();
 
-				totalQuantity.incrementValue(c, d.getProduction());
-				totalWorldInputCost.incrementValue(c, d.getProdCost());
+				if (countryNetImports > 0)
+					totalImportCommodities.incrementValue(c, countryNetImports);
+				else
+					totalExportCommodities.incrementValue(c, -countryNetImports);
 			}
-
 		}
-
-		prevWorldInputCost = totalWorldInputCost.divideBy(totalQuantity);
-		
-		// after first timestep should look at trade balance and adjust appropriately
-	/*	for (CropType crop : CropType.values()) {
-			CommodityData cd = result.getCommoditiesData().get(crop);
-			double netImport = cd.getNetImports();
-
-			if (netImport > 0)
-				totalImportCommodities.incrementValue(crop, netImport);
-			else
-				totalExportCommodities.incrementValue(crop, -netImport);
-		} */
-		
 		
+		// Look at trade balance and adjust appropriately
+		for (CropType crop : CropType.getImportedTypes()) {
+			double imports = totalImportCommodities.get(crop);
+			double exports = totalExportCommodities.get(crop);
+			double prevPrice = prevWorldPrices.get(crop);
+			double adjustedPrice = updateMarketPrices(prevPrice, imports, exports);
+			LogWriter.println(String.format("Price for %s updated from %.3f (imports %.0f, exports %.0f) to %.3f ", crop.getGamsName(), prevPrice, imports, exports, adjustedPrice));
+			prevWorldPrices.put(crop, adjustedPrice);
+		}
 		
 		// output results
 		outputTimestepResults(timestep, globalIntensityRaster, globalCropAreaRaster, globalLocationIdRaster, yieldSurfaces);
 		writeMarketFile(timestep, globalCropAreaRaster);
 	}
+	
+	public double updateMarketPrices(double previousPrice, double demand, double supply) {
+
+		if (demand > 0 || supply > 0) {
+			double ratio;
+			
+			if (demand > supply)
+				ratio = (demand-supply)/demand;
+			else
+				ratio = (demand-supply)/supply;
 
+			double adjustment = Math.exp(ratio * ModelConfig.MARKET_LAMBA);
+			return previousPrice * adjustment;
+		}
+		else {
+			return previousPrice;
+		}
+	}
 
 	private BufferedWriter getFileWriter(Timestep timestep, String file, String columnHeadings) throws IOException {
 		if (timestep.isInitialTimestep()) {
@@ -188,7 +198,7 @@ public class ModelMain {
 					timestep.getYear(), getTotalArea(LandCoverType.CROPLAND, cropAreaRaster), getTotalArea(LandCoverType.PASTURE, cropAreaRaster), getTotalArea(LandCoverType.OTHER_NATURAL, cropAreaRaster)));
 			
 			for (CropType crop : CropType.getImportedTypes() )
-				sbData.append(String.format(",%.3f", prevWorldInputCost.get(crop)));
+				sbData.append(String.format(",%.3f", prevWorldPrices.get(crop)));
 			
 			outputFile.write(sbData.toString());
 			outputFile.newLine();
@@ -314,7 +324,7 @@ public class ModelMain {
 	//			continue;
 	//		}
 
-			if (demandManager.getPopulation(country, 2010) < 8 || countryExclusionList.contains(country.getCountryName())) {
+			if (demandManager.getPopulation(country, 2010) < 80 || countryExclusionList.contains(country.getCountryName())) {
 				LogWriter.printlnError("Skipping " + country);
 				continue;
 			}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index ce2bfa74..0dc9ddf4 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -62,7 +62,7 @@ public class CountryAgent {
 		return country;
 	}
 
-	public GamsRasterOutput determineProduction(Timestep timestep, YieldRaster countryYieldSurfaces, Map<CropType, Double> worldInputEnergy) {
+	public GamsRasterOutput determineProduction(Timestep timestep, YieldRaster countryYieldSurfaces, Map<CropType, Double> worldInputPrices) {
 		currentTimestep = timestep;
 		this.countryYieldSurfaces = countryYieldSurfaces;
 		
@@ -77,7 +77,7 @@ public class CountryAgent {
 		}
 		else {
 			// optimize areas and intensity 
-			GamsRasterInput input = getGamsRasterInput(projectedDemand, worldInputEnergy);
+			GamsRasterInput input = getGamsRasterInput(projectedDemand, worldInputPrices);
 			GamsRasterOptimiser opti = new GamsRasterOptimiser(input);
 			LogWriter.println("Running " + country.getCountryName() + ", currentTimestep " + currentTimestep);
 			
@@ -89,7 +89,7 @@ public class CountryAgent {
 		return null;
 	}
 
-	private GamsRasterInput getGamsRasterInput(Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy) {
+	private GamsRasterInput getGamsRasterInput(Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputPrices) {
 
 		GamsRasterOutput prevOutput;
 		Map<CropType, Double> cropAdjs;
@@ -116,7 +116,7 @@ public class CountryAgent {
 			maxOfProdOrSupply.put(entry.getKey(), cropUsage.getProduction() + Math.max(netImports, 0));
 		}
 		
-		GamsCountryInput countryLevelInputs = GamsCountryInput.createInput(country, projectedDemand, worldInputEnergy, baseNetImport, maxOfProdOrSupply, cropAdjs, calibrate);	
+		GamsCountryInput countryLevelInputs = GamsCountryInput.createInput(country, projectedDemand, worldInputPrices, baseNetImport, maxOfProdOrSupply, cropAdjs, calibrate);	
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, prevOutput.getCropAreaRaster(), irrigationCostRaster, countryLevelInputs);
 
 		return input;
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 736cf89a..bccfbe43 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -12,7 +12,7 @@ public class GamsCountryInput {
 
 	private Country country; // not really required but useful for debugging
 	private Map<CommodityType, Double> projectedDemand;
-	private Map<CropType, Double> worldInputEnergy;
+	private Map<CropType, Double> worldInputPrices;
 	private Map<CropType, Double> maxNetImport;
 	private Map<CropType, Double> minNetImport;
 	private Map<CropType, Double> cropAdjustments;
@@ -27,19 +27,19 @@ public class GamsCountryInput {
 	private double landChangeEnergy;*/
 	
 	
-	private GamsCountryInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy,
+	private GamsCountryInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputPrices,
 			Map<CropType, Double> maxNetImport, Map<CropType, Double> minNetImport, Map<CropType, Double> cropAdjustments, boolean calibrateToObserved) {
 		super();
 		this.country = country;
 		this.projectedDemand = projectedDemand;
-		this.worldInputEnergy = worldInputEnergy;
+		this.worldInputPrices = worldInputPrices;
 		this.maxNetImport = maxNetImport;
 		this.minNetImport = minNetImport;
 		this.cropAdjustments = cropAdjustments;
 		this.calibrateToObserved = calibrateToObserved;
 	}
 	
-	public static GamsCountryInput createInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy,
+	public static GamsCountryInput createInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputPrices,
 			Map<CropType, Double> baseNetImport, Map<CropType, Double> maxOfProdOrSupply, Map<CropType, Double> cropAdjustments, boolean calibrateToObserved) {
 	
 		Map<CropType, Double> maxNetImport;
@@ -68,7 +68,7 @@ public class GamsCountryInput {
 			}
 		}
 		
-		return new GamsCountryInput(country, projectedDemand, worldInputEnergy, maxNetImport, minNetImport, cropAdjustments, calibrateToObserved);
+		return new GamsCountryInput(country, projectedDemand, worldInputPrices, maxNetImport, minNetImport, cropAdjustments, calibrateToObserved);
 	}
 	
 
@@ -80,8 +80,8 @@ public class GamsCountryInput {
 		return projectedDemand;
 	}
 
-	public Map<CropType, Double> getWorldInputEnergy() {
-		return worldInputEnergy;
+	public Map<CropType, Double> getWorldInputPrices() {
+		return worldInputPrices;
 	}
 
 	public Map<CropType, Double> getMaxNetImport() {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index d479b41c..ce2060a8 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -160,8 +160,8 @@ public class GamsLocationOptimiser {
 			}
 		}
 
-		if (DEBUG) LogWriter.println("\nWorld input energy");
-		addItemMapParm(inDB.addParameter("worldInputEnergy", 1), countryInput.getWorldInputEnergy());
+		if (DEBUG) LogWriter.println("\nWorld input prices");
+		addItemMapParm(inDB.addParameter("worldInputPrices", 1), countryInput.getWorldInputPrices());
 
 		if (DEBUG) LogWriter.println("\nMax Net Import");
 		addItemMapParm(inDB.addParameter("maxNetImport", 1), countryInput.getMaxNetImport());
-- 
GitLab