Fixes, including handling pasture and meat separately

GAMS/IntExtOpt.gms

@@ -6,6 +6,8 @@
  SET non_cereal_crop(crop_less_pasture)                                         / oilcrops, soybean, pulses, starchyRoots /;
  SET feed_crop(crop)                       / wheat, maize,                        oilcrops, soybean /;
  SET not_feed_crop(crop)                                 / rice, tropicalCereals,                    pulses, starchyRoots, pasture /;
+ SET import_crop(all_types) / meat,          wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots /;
+ 
  
  SET location;
  PARAMETER suitableLandArea(location)        areas of land in Mha;
@@ -17,9 +19,9 @@
  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(crop)            average input energy from world exports used to determine if we should import or export;
- PARAMETER maxNetImport(crop)                maximum net import for each crop based on world market;
- PARAMETER minNetImport(crop)                minimum net import for each crop based on world market;
+ PARAMETER worldInputEnergy(all_types)       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;
  
  SCALAR meatEfficency                efficiency of converting feed and pasture into animal products;
@@ -35,17 +37,15 @@ $load meatEfficency, maxLandUseChange, minFeedRate, irrigCost
 $gdxin
   
  SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /  
- 
-* display demand;
-    
+     
  SCALAR minDemandPerCereal / 0.1 /
-               
- PARAMETER feedDM(crop) energy from feed in MJ per t (wet) (need to check these figures)
-          /   wheat            0.87
-              maize            0.87
-              oilcrops         0.32
-              soybean          0.87
-              pasture          0.2 / ;
+         
+ PARAMETER feedDM(crop) energy from feed in MJ per kg (wet) (i.e. MJ per kg DM multiplied by DM kg per kg)
+          /   wheat     11.7   
+              maize     11.1  
+              oilcrops  10.8  
+              soybean   10.6  
+              pasture   5.3    / ; 
  
  VARIABLES
        area(crop, location)               total area for each crop - Mha
@@ -53,7 +53,7 @@ $gdxin
        irrigI(crop, location)             irrigation intensity for each crop - factor between 0 and 1
        otherIntensity(crop, location)     other intensity for each crop - unit less
        feedAmount(crop)                   amount of feed for each crop - Mt
-       netImportAmount(all_types)         net imports of crops and meat - Mt
+       netImportAmount(all_types)       net imports of crops and meat - Mt
        yield(crop, location)              yield per area for each crop - t per ha
        unitEnergy(crop, location)         energy per area for each crop - energy
        net_supply(crop_less_pasture)      supply after exports and feed
@@ -75,21 +75,22 @@ $gdxin
        MEAT_DEMAND_CONSTRAINT                           satisfy demand for meat
        MAX_FERT_INTENSITY_CONSTRAINT(crop, location)    constraint on maximum fertilizer intensity
        MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location)   constraint on maximum irrigation intensity
+*       MAX_OTHER_INTENSITY_CONSTRAINT(crop, location)   constraint on maximum other intensity
        TOTAL_LAND_CHANGE_CONSTRAINT(location)           constraint on the rate of land use change 
        CROPLAND_CHANGE_CONSTRAINT(location)             constraint on the rate of land use change 
        PASTURE_CHANGE_CONSTRAINT(location)              constraint on the rate of land use change 
        AGRI_LAND_INCREASE_CONSTRAINT(location)          constraint expansion of agricultural land
        AGRI_LAND_DECREASE_CONSTRAINT(location)          constraint expansion of agricultural land
        NON_FEED_CROP_CONSTRAINT(not_feed_crop)          constraint to set non feed crop feed usage to zero
-       MAX_NET_IMPORT_CONSTRAINT(crop)                  constraint on max net imports
-       MIN_NET_IMPORT_CONSTRAINT(crop)                  constraint on min net imports
+       MAX_NET_IMPORT_CONSTRAINT(import_crop)           constraint on max net imports
+       MIN_NET_IMPORT_CONSTRAINT(import_crop)           constraint on min net imports
        MIN_FEED_CONSTRAINT                              constraint on min feed rate
        NET_SUPPLY_EQ(crop_less_pasture)                 calc net supply for crops
        ENERGY_EQ                                        total energy objective function;
  
  UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= fertI(crop, location) + 
-                                                                       (1+irrigI(crop, location)*irrigCost(location) ** 1.3) +
-                                                                       (1+otherIntensity(crop, location) ** 2);
+                                                                     ((1+irrigI(crop, location) ** 2) * irrigCost(location)) +
+                                                                     (1+otherIntensity(crop, location) ** 2);
  
  YIELD_EQ(crop, location) .. yield(crop, location) =E= (
                yieldNone(crop, location) + 
@@ -119,6 +120,7 @@ $gdxin
     
  MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
  MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
+* MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) .. otherIntensity(crop, location) =L= 10;
  
  TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop, area(crop, location));
  
@@ -132,20 +134,20 @@ $gdxin
  
  NON_FEED_CROP_CONSTRAINT(not_feed_crop) .. feedAmount(not_feed_crop) =E= 0;
  
- MAX_NET_IMPORT_CONSTRAINT(crop) .. netImportAmount(crop) =L= maxNetImport(crop);
+ MAX_NET_IMPORT_CONSTRAINT(import_crop) .. netImportAmount(import_crop) =L= maxNetImport(import_crop);
  
- MIN_NET_IMPORT_CONSTRAINT(crop) .. netImportAmount(crop) =G= minNetImport(crop);
+ MIN_NET_IMPORT_CONSTRAINT(import_crop) .. netImportAmount(import_crop) =G= minNetImport(import_crop);
  
  MIN_FEED_CONSTRAINT .. sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= minFeedRate * (demand('meat') - netImportAmount('meat'));
  
  ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location))
             + sum(location, sum(crop, area(crop, location)) - sum(crop, previousArea(crop, location)) * landChangeEnergy)
-            + sum(crop, (netImportAmount(crop)) * worldInputEnergy(crop))) / 1000000;
+            + sum(import_crop, (netImportAmount(import_crop)) * worldInputEnergy(import_crop))) / 1000000;
  
  MODEL LAND_USE /ALL/ ;
  SOLVE LAND_USE USING NLP MINIMIZING energy;
    
-display net_supply.l, area.l, yield.l, feedAmount.l, netImportAmount.l
+display feedAmount.l, minFeedRate, demand, netImportAmount.l;
     
  Scalar ms 'model status', ss 'solve status'; 
  ms=LAND_USE.modelstat;

src/ac/ed/lurg/ModelConfig.java

@@ -93,7 +93,7 @@ public class ModelConfig {
 	public static final String BASELINE_CONSUMP_FILE = DATA_DIR + File.separator + "base_consump.csv";
 	public static final String COUNTRY_CODES_FILE = DATA_DIR + File.separator + "country_codes3.csv";
 	public static final String COUNTRY_DATA_FILE = DATA_DIR + File.separator + "country_data.csv";
-	public static final String COMMODITY_DATA_FILE = DATA_DIR + File.separator + "con_prod_c.csv";
+	public static final String COMMODITY_DATA_FILE = DATA_DIR + File.separator + "con_prod_c_and_m.csv";
 
 	// yield data
 	public static final String YIELD_DIR = "/Users/peteralexander/Documents/LURG/LPJ/tom-April2015";

src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java

@@ -60,6 +60,9 @@ public class GamsLocationOptimiser {
 
 		GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.GAMS_MODEL);
 		GAMSOptions opt = ws.addOptions();
+		
+		opt.setAllModelTypes("conopt");
+		
 		opt.defines("gdxincname", inDB.getName());
 
 		long startTime = System.currentTimeMillis();
@@ -128,7 +131,7 @@ public class GamsLocationOptimiser {
 			String locString = Integer.toString(locationId);
 			YieldResponsesItem yresp = entry.getValue();
 
-			for (CropType crop : CropType.values()) {
+			for (CropType crop : CropType.getNonMeatTypes()) {
 				if (DEBUG) LogWriter.println(String.format("     %15s,\t %.1f,\t %.1f, \t %.1f,\t %.1f,\t\t [%.2f],\t [%.2f]", 
 						crop.getGamsName(), yresp.getYieldNone(crop), yresp.getYieldFertOnly(crop), yresp.getYieldIrrigOnly(crop), yresp.getYieldMax(crop), yresp.getFertParam(crop), yresp.getIrrigParam(crop)));
 
@@ -205,8 +208,8 @@ public class GamsLocationOptimiser {
 
 			if (!cropType.equals(prevCropType)) {
 				feedAmount = varFeedAmount.findRecord(itemName).getLevel();
-				netImport = varNetImports.findRecord(itemName).getLevel();
-
+				netImport = cropType.isImportedCrop() ? varNetImports.findRecord(itemName).getLevel() : 0;
+					
 				cropUsageData.put(cropType, new CropUsageData(feedAmount, netImport));
 				if (DEBUG) LogWriter.println(String.format("\n%s:\tfeedAmount= %.1f,\tnetImports= %.3f", itemName, feedAmount, netImport)); 
 			}

src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java

@@ -191,8 +191,8 @@ public class GamsRasterOptimiser {
 
 		}
 
-		int numCerealCats = 1;
-		int numPastureCats = 1;
+		int numCerealCats = 2;
+		int numPastureCats = 2;
 
 		int thisShouldLookAtCropsOtherThanJustWheat; // need to consider other crops, and perhaps other yieldTypes as well
 		List<Double> wheatlDivisions = getDivisions(yieldRaster, CropType.WHEAT, numCerealCats); 
@@ -253,7 +253,7 @@ public class GamsRasterOptimiser {
 				}
 
 				// Crop yields
-				for (CropType crop : CropType.values()) {
+				for (CropType crop : CropType.getNonMeatTypes()) {
 					for (YieldType yieldType : YieldType.values()) {
 						avgYResp.setYield(yieldType, crop, yresp.getYield(yieldType, crop));
 					}

src/ac/ed/lurg/types/CropType.java

@@ -2,6 +2,7 @@ package ac.ed.lurg.types;
 
 import java.util.Collection;
 import java.util.HashMap;
+import java.util.HashSet;
 import java.util.Map;
 
 import ac.ed.lurg.utils.LogWriter;
@@ -16,16 +17,45 @@ public enum CropType {
 	SOYBEAN("Soyabeans", "soybean"),
 	PULSES("Pulses + (Total)", "pulses"),
 	STARCHY_ROOTS("Starchy Roots + (Total)", "starchyRoots"),
-	PASTURE("meatmilkeggs", "pasture");
+	MEAT("meatmilkeggs", "meat", false, true),
+	PASTURE("pasture", "pasture", false, false);
 
 	private String faoName;
 	private String gamsName;
-	
-	CropType (String faoName, String gamsName) {
+	private boolean importedCrop;
+	private boolean isMeat;
+
+	CropType (String faoName, String gamsName, boolean importedCrop, boolean isMeat) {
 		this.faoName = faoName;
 		this.gamsName = gamsName;
+		this.importedCrop = importedCrop;
+		this.isMeat = isMeat;
+	}
+	
+	CropType (String faoName, String gamsName) {
+		this(faoName, gamsName, true, false);
+	}
+	
+	public static Collection<CropType> getImportedTypes() {
+		Collection<CropType> comms = new HashSet<CropType>();
+		
+		for (CropType c : values())
+			if (c.importedCrop)
+				comms.add(c);
+		
+		return comms;
 	}
 	
+	public static Collection<CropType> getNonMeatTypes() {
+		Collection<CropType> comms = new HashSet<CropType>();
+		
+		for (CropType c : values())
+			if (!c.isMeat)
+				comms.add(c);
+		
+		return comms;
+	}
+
 	private static final Map<String, CropType> faoNameCache = new HashMap<String, CropType>();
 	private static final Map<String, CropType> gamsNameCache = new HashMap<String, CropType>();
     static {
@@ -67,4 +97,8 @@ public enum CropType {
 	public String getGamsName() {
 		return gamsName;
 	}
+	
+	public boolean isImportedCrop() {
+		return importedCrop;
+	}
 }

src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java

@@ -55,28 +55,29 @@ public class LPJYieldResponseMapReader {
 						record.y = d;
 						break;
 					case 2:
-						record.teSW = d;
+						//YieldRecord values in kg DM / m2, we want t/ha, so 10 times larger
+						record.teSW = d * 10;
 						break;
 					case 3:
-						record.teSWirr = d;
+						record.teSWirr = d * 10;
 						break;
 					case 4:
-						record.teWW = d;
+						record.teWW = d * 10;
 						break;
 					case 5:
-						record.teWWirr = d;
+						record.teWWirr = d * 10;
 						break;
 					case 6:
-						record.teCo = d;
+						record.teCo = d * 10;
 						break;
 					case 7:
-						record.teCoirr = d;
+						record.teCoirr = d * 10;
 						break;
 					case 8:
-						record.trRi = d;
+						record.trRi = d * 10;
 						break;
 					case 9:
-						record.trRiirr = d;
+						record.trRiirr = d * 10;
 						break;
 					}
 
@@ -105,11 +106,25 @@ public class LPJYieldResponseMapReader {
 	
 	public void setData(YieldType noIrrigYieldType, YieldType maxIrrigYieldType, YieldResponsesItem item, YieldRecord record) {
 		
-		//YieldRecord values in kg DM / m2, we want t/ha, so 10 times larger
-		for (CropType crop : CropType.getAllItems()) {
-			item.setYield(noIrrigYieldType, crop, record.teSW * 10);
-			item.setYield(maxIrrigYieldType, crop, record.teSWirr * 10);
-		}
+		item.setYield(noIrrigYieldType, CropType.WHEAT, Math.max(record.teSW, record.teWW));
+		item.setYield(maxIrrigYieldType, CropType.WHEAT, Math.max(record.teSWirr, record.teWWirr));
+		item.setYield(noIrrigYieldType, CropType.MAIZE, record.teCo);
+		item.setYield(maxIrrigYieldType, CropType.MAIZE, record.teCoirr);
+		item.setYield(noIrrigYieldType, CropType.RICE, record.trRi);
+		item.setYield(maxIrrigYieldType, CropType.RICE, record.trRiirr);
+		
+		item.setYield(noIrrigYieldType, CropType.TROPICAL_CEREALS, record.teCo);
+		item.setYield(maxIrrigYieldType, CropType.TROPICAL_CEREALS, record.teCoirr);
+		item.setYield(noIrrigYieldType, CropType.OILCROPS, record.teWW);
+		item.setYield(maxIrrigYieldType, CropType.OILCROPS, record.teWWirr);
+		item.setYield(noIrrigYieldType, CropType.SOYBEAN, record.teSW);
+		item.setYield(maxIrrigYieldType, CropType.SOYBEAN, record.teSWirr);
+		item.setYield(noIrrigYieldType, CropType.PULSES, record.teSW);
+		item.setYield(maxIrrigYieldType, CropType.PULSES, record.teSWirr);
+		item.setYield(noIrrigYieldType, CropType.STARCHY_ROOTS, record.teSW * 1.1);
+		item.setYield(maxIrrigYieldType, CropType.STARCHY_ROOTS, record.teSWirr * 1.1);
+		item.setYield(noIrrigYieldType, CropType.PASTURE, record.teSW * 1.2);
+		item.setYield(maxIrrigYieldType, CropType.PASTURE, record.teSWirr * 1.2);
 	}
 	
 	class YieldRecord {

src/ac/ed/lurg/yield/YieldResponse.java

@@ -4,7 +4,6 @@ import java.util.HashMap;
 import java.util.Map;
 
 import ac.ed.lurg.types.YieldType;
-import ac.ed.lurg.utils.LogWriter;
 
 public class YieldResponse {
 	Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
@@ -21,7 +20,7 @@ public class YieldResponse {
 		Double d = yields.get(type);
 		
 		if (d == null) {
-			LogWriter.printlnError("YieldResponse: No yield found for type " + type);
+		//	LogWriter.printlnError("YieldResponse: No yield found for type " + type);
 			return Double.NaN;
 		}
 		
@@ -30,14 +29,14 @@ public class YieldResponse {
 	
 	public double getFertParam() {
 		if (fertParm == 0)
-			fertParm = calcParam (0, 0.6, 1); // we do have MID fert data, but looks wrong calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG));
+			fertParm = calcParam (0, 0.7, 1); // we do have MID fert data, but looks wrong calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG));
 		
 		return fertParm;
 	}
 	
 	public double getIrrigParam() {
 		if (irrigParm == 0)
-			irrigParm = calcParam (0, 0.6, 1);  // we don't have a mid irrigation figure, so lets assume 60% at mid point
+			irrigParm = calcParam (0, 0.7, 1);  // we don't have a mid irrigation figure, so lets assume 60% at mid point
 		
 		return irrigParm;
 	}

src/ac/ed/lurg/yield/YieldResponsesItem.java

@@ -5,6 +5,7 @@ import java.util.Map;
 
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.YieldType;
+import ac.ed.lurg.utils.LogWriter;
 import ac.sac.raster.RasterItem;
 
 public class YieldResponsesItem implements RasterItem {
@@ -26,6 +27,9 @@ public class YieldResponsesItem implements RasterItem {
 	public double getYield(YieldType yieldType, CropType crop) {
 		YieldResponse yr = getYieldResponseForCrop(crop);
 		double d = yr.getYield(yieldType);
+		if (Double.isNaN(d)) {
+			LogWriter.printlnError("YieldResponses: No yield found for type " + yieldType + ", crop " + crop);
+		}
 		return d;
 	}