From eab29472888bcda03ea540760901b4890133096b Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Sun, 3 Jul 2016 22:01:56 +0100
Subject: [PATCH] Some simplifications

---
 GAMS/IntExtOpt.gms              | 51 ++++++++++++++++-----------------
 src/ac/ed/lurg/ModelConfig.java | 10 +++----
 src/ac/ed/lurg/ModelMain.java   | 11 +++++--
 3 files changed, 38 insertions(+), 34 deletions(-)

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 26fb1bbf..fe074b39 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -3,10 +3,10 @@
  SET crop(all_types)                       / wheat, maize, rice, oilcrops, pulses, starchyRoots, pasture /;
  SET crop_less_pasture(crop)               / wheat, maize, rice, oilcrops, pulses, starchyRoots /;
  SET cereal_crop(crop)                     / wheat, maize, rice /;
- SET non_cereal_crop(crop)                                     / oilcrops, pulses, starchyRoots,              pasture /;
+ SET non_cereal_crop(crop)                                     / oilcrops, pulses, starchyRoots, pasture /;
  SET feed_crop(crop)                       / wheat, maize,       oilcrops,                       pasture/;
  SET feed_crop_less_pasture(feed_crop)     / wheat, maize,       oilcrops          /;
- SET not_feed_crop(crop)                                 / rice,                    pulses, starchyRoots /;
+ SET not_feed_crop(crop)                                 / rice,           pulses, starchyRoots /;
  SET import_crop(all_types) / meat,          wheat, maize, rice, oilcrops, pulses, starchyRoots /;
  
  SET location;
@@ -48,25 +48,24 @@ $gdxin
  SCALAR unhandledCropArea includes fruit veg forage crops set aside and failed crop / 0.4 /;
  previousArea(crop_less_pasture, location) = previousArea(crop_less_pasture, location) * (1.0 - unhandledCropArea);
  
- PARAMETER feedDM(crop)  kg DM per kg of feed the conversion from feed to meat is done in the R animal product index
-          /   wheat     0.89   
-              maize     0.7  
-              oilcrops  0.9  
-              pasture   1    / ; 
-* pasture yield is done in DM terms
+ PARAMETER cropDM(crop)  kg DM per kg of feed the conversion from feed to meat is done in the R animal product index
+          /   wheat         0.87   
+              maize         0.86  
+              rice          0.89
+              oilcrops      0.88  
+              pulses        0.31
+              starchyRoots  0.21
+              pasture       1    / ; 
    
  PARAMETER baseCost(crop)  cost per ha before intensity values not including fertiliser or irrigation i.e. seed spray and other.  In 1000 $ per ha
-          /   wheat             1.1   
-              maize             0.9 
-              rice              1.1 
-              oilcrops          0.7 
-              pulses            0.5
-              starchyRoots      4.0
+          /   wheat             0.31   
+              maize             0.26 
+              rice              0.4
+              oilcrops          0.24
+              pulses            0.31
+              starchyRoots      4.64
               pasture           0.0  / ; 
-               
-* wheat             1.4 * 0.223   
-*              oilcrops          1.4 * 0.169 
-               
+                              
  VARIABLES
        area(crop, location)               total area for each crop - Mha
        fertI(crop, location)              fertilizer intensity for each crop - factor between 0 and 1
@@ -117,7 +116,7 @@ $gdxin
  UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E=  ( baseCost(crop) +              
                                                                      fertiliserUnitCost * fertI(crop, location) + 
                                                                      irrigCost(location) * irrigMaxRate(crop, location) * irrigI(crop, location) +
-                                                                     (otherIntensity(crop, location))
+                                                                     otherIntensity(crop, location)
                                                                      ) ;
  
  YIELD_EQ(crop, location) .. yield(crop, location) =E= (
@@ -132,17 +131,17 @@ $gdxin
                (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
                                       (1 - exp(-fertParam(crop, location)*fertI(crop, location))) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location)))
 
-            ) * cropAdj(crop) * (1 - exp(-otherIntensity(crop, location)*5));
+            ) * cropAdj(crop) * (1 - exp(-otherIntensity(crop, location)*4));
    
- NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + importAmount(crop) - exportAmount(crop);
+ NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= (sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop)) / cropDM(crop) + importAmount(crop) - exportAmount(crop);
   
- NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop) ;
+ NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop);
   
  CEREAL_DEMAND_CONSTRAINT(cereal_crop) .. net_supply(cereal_crop) =G= demand('cereals') * minDemandPerCereal;
   
  TOTAL_CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, net_supply(cereal_crop)) =G= demand('cereals');
  
- MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= demand('meat') - importAmount('meat') + exportAmount('meat');
+ MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedAmount(feed_crop)) =G= demand('meat') - importAmount('meat') + exportAmount('meat');
     
  MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
  MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
@@ -157,7 +156,7 @@ $gdxin
  PASTURE_IMPORT_CONSTRAINT .. importAmount('pasture') =E= 0;
  PASTURE_EXPORT_CONSTRAINT ..  exportAmount('pasture') =E= 0;
  
- MIN_FEED_CONSTRAINT .. sum(feed_crop_less_pasture, feedDM(feed_crop_less_pasture) * feedAmount(feed_crop_less_pasture)) =G= minFeedRate * (demand('meat') - importAmount('meat') + exportAmount('meat'));
+ MIN_FEED_CONSTRAINT .. sum(feed_crop_less_pasture, feedAmount(feed_crop_less_pasture)) =G= minFeedRate * (demand('meat') - importAmount('meat') + exportAmount('meat'));
  
  IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * area(crop, location)) / suitableLandArea(location);
  
@@ -186,7 +185,7 @@ $gdxin
    
  fertI.L(crop, location) = 0.5;
  irrigI.L(crop, location) = 0.5;
- otherIntensity.L(crop, location) = 0.1;
+ otherIntensity.L(crop, location) = 0.5;
 
  importAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) gt 0) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
  exportAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) lt 0) = -(maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
@@ -207,7 +206,7 @@ $gdxin
  parameter feedEnergy(all_types);
  
  totalProd(crop) = sum(location, area.l(crop, location) * yield.l(crop, location));
- totalProd('meat') = meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount.l(feed_crop));
+ totalProd('meat') = meatEfficency*sum(feed_crop, feedAmount.l(feed_crop));
  totalProdCost(crop) = sum(location, unitEnergy.l(crop, location) * area.l(crop, location));
  totalCropland(location) = sum(crop_less_pasture, area.l(crop_less_pasture, location));
  
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index c0a66513..7bd58e23 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -139,17 +139,17 @@ public class ModelConfig {
 	public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", 0.1);
 
 	public static final double PASTURE_HARVEST_FRACTION = getDoubleProperty("PASTURE_HARVEST_FRACTION", 0.25);
-	public static final double MEAT_EFFICIENCY = getDoubleProperty("MEAT_EFFICIENCY", 1.0);
-	public static final double IRRIGIATION_EFFICIENCY = getDoubleProperty("IRRIGIATION_EFFICIENCY", 0.4);
-	public static final double LAND_CHANGE_COST = getDoubleProperty("LAND_CHANGE_COST", 2.0);
+	public static final double MEAT_EFFICIENCY = getDoubleProperty("MEAT_EFFICIENCY", 1.0);  // 'meat' is includes feed conversion ratio already, this is tech. change or similar
+	public static final double IRRIGIATION_EFFICIENCY = getDoubleProperty("IRRIGIATION_EFFICIENCY", 0.5);
+	public static final double LAND_CHANGE_COST = getDoubleProperty("LAND_CHANGE_COST", 1.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.15);  
 
 	public static final double MARKET_LAMBA = getDoubleProperty("MARKET_LAMBA", 0.3); // controls international market price adjustment rate
 	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", 0.02);
-	public static final double FERTILISER_COST_PER_T = getDoubleProperty("FERTILISER_COST_PER_T", 5.0);
+	public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 0.01);
+	public static final double FERTILISER_COST_PER_T = getDoubleProperty("FERTILISER_COST_PER_T", 0.87 * 1.4 * 2); // £900/t N *  1.4$/£ * 2NPK/N
 	public static final double FERTILISER_MAX_COST = FERTILISER_COST_PER_T * MAX_FERT_AMOUNT/1000;
 	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.3);  // price factor in international trade, transport cost and real trade barriers
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 79a51abc..beaf1e07 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -78,8 +78,13 @@ public class ModelMain {
 
 		// in first timestep we don't have this info, but ok as constrained to import/export specified amount
 		prevWorldPrices = new HashMap<CropType, GlobalPrice>();  
-		for (CropType c : CropType.getImportedTypes()) 
-			prevWorldPrices.put(c, GlobalPrice.createInitial(0.5));
+		prevWorldPrices.put(CropType.WHEAT, GlobalPrice.createInitial(0.15));
+		prevWorldPrices.put(CropType.MAIZE, GlobalPrice.createInitial(0.15));
+		prevWorldPrices.put(CropType.RICE, GlobalPrice.createInitial(0.4));
+		prevWorldPrices.put(CropType.OILCROPS, GlobalPrice.createInitial(0.4));
+		prevWorldPrices.put(CropType.PULSES, GlobalPrice.createInitial(0.2));
+		prevWorldPrices.put(CropType.STARCHY_ROOTS, GlobalPrice.createInitial(0.1));
+		prevWorldPrices.put(CropType.MEAT, GlobalPrice.createInitial(0.2));
 	}
 
 	/* run the model */
@@ -336,7 +341,7 @@ public class ModelMain {
 
 			// DEBUG code
 			if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {
-				if (!(cc.getName().equals("United States of America") || cc.getName().equals("Chinaxx") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
+				if (!(cc.getName().equals("United States of Americaxxx") || cc.getName().equals("Chinaxx") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("Philippines")) ) {
 					continue;
 				}
 			}
-- 
GitLab