Split pasture and meat and fixes

GAMS/IntExtOpt.gms

- SET all_types      / cereals, wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots, pasture_or_meat /;
+ SET all_types              / meat, cereals, wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots, pasture /;
  
- SET crop(all_types)         / wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots, pasture_or_meat /;
- SET crop_less_pasture(crop) / wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots /;
- SET cereal_crop(crop)       / wheat, maize, rice, tropicalCereals /;
- SET non_cereal_crop(crop)                                        / oilcrops, soybean, pulses, starchyRoots /;
- SET feed_crop(crop)         / wheat, maize,                        oilcrops, soybean /;
- SET not_feed_crops(crop)                  / rice, tropicalCereals,                    pulses, starchyRoots, pasture_or_meat /;
+ SET crop(all_types)                       / wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots, pasture /;
+ SET crop_less_pasture(crop)               / wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots /;
+ SET cereal_crop(crop_less_pasture)        / wheat, maize, rice, tropicalCereals /;
+ 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 location;
  PARAMETER suitableLandArea(location)        areas of land in Mha;
@@ -25,7 +25,7 @@
  SCALAR meatEfficency                efficiency of converting feed and pasture into animal products;
  SCALAR maxLandUseChange             max rate of land use change;
  SCALAR landChangeEnergy             energy required to add ha of agricultural land;
- SCALAR minFeedRate                  minimum rate of feed for producing animal products;
+ SCALAR minFeedRate                  minimum rate of feed for producing animal products (why is this needed?);
 
 $gdxin %gdxincname%
 $load location, suitableLandArea, previousArea, demand, landChangeEnergy
@@ -36,14 +36,16 @@ $gdxin
   
  SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /  
  
- display previousArea
+* 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_or_meat  1.0 / ;
+              pasture          0.2 / ;
  
  VARIABLES
        area(crop, location)               total area for each crop - Mha
@@ -51,9 +53,10 @@ $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(crop)              export of crop - 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
        energy                             total input energy - energy;
  
  POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount;
@@ -67,7 +70,8 @@ $gdxin
        UNIT_ENERGY_EQ(crop, location)                   energy per area
        YIELD_EQ(crop, location)                         yield given chosen intensity
        NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop)    satisfy demand for non-cereal crop 
-       CEREAL_DEMAND_CONSTRAINT                         satisfy demand for cereal
+       CEREAL_DEMAND_CONSTRAINT(cereal_crop)            satisfy demand for cereal so that exports can at least be met. Could also allow min. proporation of cereal consumption of each type
+       TOTAL_CEREAL_DEMAND_CONSTRAINT                   satisfy demand for cereal
        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
@@ -76,10 +80,11 @@ $gdxin
        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_crops)         constraint to set non feed crop feed usage to zero
+       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
        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) + 
@@ -101,36 +106,37 @@ $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)))
 *            ) * otherIntensity(crop, location);
   
- NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. sum(location, area(non_cereal_crop, location) * yield(non_cereal_crop, location)) - feedAmount(non_cereal_crop) =G= 
-            demand(non_cereal_crop) - netImportAmount(non_cereal_crop);
- 
- CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, sum(location, area(cereal_crop, location) * yield(cereal_crop, location)) - feedAmount(cereal_crop)) =G= 
-            demand('cereals') - sum(cereal_crop, netImportAmount(cereal_crop));
+ NET_SUPPLY_EQ(crop_less_pasture) .. net_supply(crop_less_pasture) =E= sum(location, area(crop_less_pasture, location) * yield(crop_less_pasture, location)) - feedAmount(crop_less_pasture) + netImportAmount(crop_less_pasture);
+  
+ 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)) + feedDM('pasture_or_meat') * sum(location, area('pasture_or_meat', location) * yield('pasture_or_meat', location)) =G= 
-            demand('pasture_or_meat') - netImportAmount('pasture_or_meat');
+ MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) + feedDM('pasture') * sum(location, area('pasture', location) * yield('pasture', location)) =G= 
+            demand('meat') - netImportAmount('meat');
     
  MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
- 
  MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
  
  TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop, area(crop, location));
  
  CROPLAND_CHANGE_CONSTRAINT(location) .. sum(crop_less_pasture, area(crop_less_pasture, location)) =L= (1 + maxLandUseChange) * sum(crop_less_pasture, previousArea(crop_less_pasture, location));
  
- PASTURE_CHANGE_CONSTRAINT(location) ..  area('pasture_or_meat', location) =L= (1 + maxLandUseChange) * previousArea('pasture_or_meat', location);
+ PASTURE_CHANGE_CONSTRAINT(location) ..  area('pasture', location) =L= (1 + maxLandUseChange) * previousArea('pasture', location);
  
  AGRI_LAND_INCREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) =L= (1 + maxLandUseChange) * sum(crop, previousArea(crop, location));
  
  AGRI_LAND_DECREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) =G= (1 - maxLandUseChange) * sum(crop, previousArea(crop, location));
  
- NON_FEED_CROP_CONSTRAINT(not_feed_crops) .. feedAmount(not_feed_crops) =E= 0;
+ NON_FEED_CROP_CONSTRAINT(not_feed_crop) .. feedAmount(not_feed_crop) =E= 0;
  
  MAX_NET_IMPORT_CONSTRAINT(crop) .. netImportAmount(crop) =L= maxNetImport(crop);
  
  MIN_NET_IMPORT_CONSTRAINT(crop) .. netImportAmount(crop) =G= minNetImport(crop);
  
- MIN_FEED_CONSTRAINT .. sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= minFeedRate * (demand('pasture_or_meat') - netImportAmount('pasture_or_meat'));
+ 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)
@@ -138,7 +144,9 @@ $gdxin
  
  MODEL LAND_USE /ALL/ ;
  SOLVE LAND_USE USING NLP MINIMIZING energy;
-  
+   
+display net_supply.l, area.l, yield.l, feedAmount.l, netImportAmount.l
+    
  Scalar ms 'model status', ss 'solve status'; 
  ms=LAND_USE.modelstat;
  ss=LAND_USE.solvestat;

src/ac/ed/lurg/country/CountryAgent.java

@@ -56,7 +56,7 @@ public class CountryAgent {
 			
 			if (landCover != null) {
 				areasItem.setCropArea(CropType.WHEAT, landCover.getLandCover(LandDataType.CROPLAND));  // we allow free substitution between crop types so this is ok-ish
-				areasItem.setCropArea(CropType.MEAT_OR_PASTURE, landCover.getLandCover(LandDataType.PASTURE));
+				areasItem.setCropArea(CropType.PASTURE, landCover.getLandCover(LandDataType.PASTURE));
 				areasItem.setLandCoverArea(LandDataType.FOREST, landCover.getLandCover(LandDataType.FOREST));
 				areasItem.setLandCoverArea(LandDataType.OTHER_NATURAL, landCover.getLandCover(LandDataType.OTHER_NATURAL));
 			}

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

@@ -254,9 +254,9 @@ public class GamsLocationOptimiser {
 	private void addCommodityMapParm(GAMSParameter parm, Map<CommodityType, Double> itemMap) {
 		for (Map.Entry<CommodityType, Double> entry : itemMap.entrySet()) {
 			double d = entry.getValue();
-			if (DEBUG) LogWriter.println(String.format("     %15s,\t %.1f", entry.getKey().getFaoName(), d));
+			if (DEBUG) LogWriter.println(String.format("     %15s,\t %.1f", entry.getKey().getGamsName(), d));
 			if (!Double.isNaN(d))
-				parm.addRecord(entry.getKey().getFaoName()).setValue(d);
+				parm.addRecord(entry.getKey().getGamsName()).setValue(d);
 		}
 	}
 

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

@@ -35,7 +35,7 @@ public class GamsLocationTest {
 		dummyMap.put(CommodityType.OILCROPS, 50.0);
 		dummyMap.put(CommodityType.PULSES, 60.0);
 		dummyMap.put(CommodityType.STARCHY_ROOTS, 150.0);
-		dummyMap.put(CommodityType.MEAT_OR_PASTURE, 480.0);
+		dummyMap.put(CommodityType.MEAT, 480.0);
 		return dummyMap;
 	}
 
@@ -56,7 +56,7 @@ public class GamsLocationTest {
 		
 		for (CropType crop : CropType.getAllItems()) {
 			double factor = base;
-			if (CropType.MEAT_OR_PASTURE == crop)
+			if (CropType.PASTURE == crop)
 				factor = base * pastureFactor;
 				
 			yresp.setYield(YieldType.NO_FERT_NO_IRRIG, crop, factor);
@@ -102,7 +102,7 @@ public class GamsLocationTest {
 		dummyMap.setCropArea(CropType.OILCROPS, 5.0);
 		dummyMap.setCropArea(CropType.PULSES, 5.0);
 		dummyMap.setCropArea(CropType.STARCHY_ROOTS, 5.0);
-		dummyMap.setCropArea(CropType.MEAT_OR_PASTURE, 5.0 + pasture);
+		dummyMap.setCropArea(CropType.PASTURE, 5.0 + pasture);
 		
 		dummyMap.setLandCoverArea(LandDataType.FOREST, 200 - pasture - cereal);
 		dummyMap.setLandCoverArea(LandDataType.OTHER_NATURAL, 30);
@@ -119,7 +119,7 @@ public class GamsLocationTest {
 		dummyMap.put(CropType.OILCROPS, 10.0);
 		dummyMap.put(CropType.PULSES, 20.0);
 		dummyMap.put(CropType.STARCHY_ROOTS, 5.0);
-		dummyMap.put(CropType.MEAT_OR_PASTURE, 80.0);
+		dummyMap.put(CropType.PASTURE, 80.0);
 		return dummyMap;
 	}
 
@@ -133,7 +133,7 @@ public class GamsLocationTest {
 		dummyMap.put(CropType.OILCROPS, 10.0);
 		dummyMap.put(CropType.PULSES, 20.0);
 		dummyMap.put(CropType.STARCHY_ROOTS, 5.0);
-		dummyMap.put(CropType.MEAT_OR_PASTURE, 5.0);
+		dummyMap.put(CropType.PASTURE, 5.0);
 		return dummyMap;
 	}
 
@@ -147,7 +147,7 @@ public class GamsLocationTest {
 		dummyMap.put(CropType.OILCROPS, -10.0);
 		dummyMap.put(CropType.PULSES, -20.0);
 		dummyMap.put(CropType.STARCHY_ROOTS, -5.0);
-		dummyMap.put(CropType.MEAT_OR_PASTURE, -5.0);
+		dummyMap.put(CropType.PASTURE, -5.0);
 		return dummyMap;
 	}
 }

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

@@ -191,12 +191,12 @@ public class GamsRasterOptimiser {
 
 		}
 
-		int numCerealCats = 3;
-		int numPastureCats = 3;
+		int numCerealCats = 1;
+		int numPastureCats = 1;
 
 		int thisShouldLookAtCropsOtherThanJustWheat; // need to consider other crops, and perhaps other yieldTypes as well
 		List<Double> wheatlDivisions = getDivisions(yieldRaster, CropType.WHEAT, numCerealCats); 
-		List<Double> pastureDivisions = getDivisions(yieldRaster, CropType.MEAT_OR_PASTURE, numPastureCats);
+		List<Double> pastureDivisions = getDivisions(yieldRaster, CropType.PASTURE, numPastureCats);
 
 		if (DEBUG) LogWriter.println("Making " + numCerealCats * numPastureCats + " categories");
 
@@ -236,7 +236,7 @@ public class GamsRasterOptimiser {
 				IrrigationCostItem irrigCost = irrigCostRaster.get(key);
 
 				int cerealCat = findCategory(wheatlDivisions, yresp.getYieldMax(CropType.WHEAT));
-				int pastureCat = findCategory(pastureDivisions, yresp.getYieldMax(CropType.MEAT_OR_PASTURE));
+				int pastureCat = findCategory(pastureDivisions, yresp.getYieldMax(CropType.PASTURE));
 				Integer id = cerealCat + pastureCat * numCerealCats;
 
 				AveragingYieldResponsesItem avgYResp = aggregatedYields.lazyGet(id);
@@ -278,7 +278,7 @@ public class GamsRasterOptimiser {
 		for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) {
 			LogWriter.println(e.getKey() + " category has " + e.getValue().size() + " raster areas");
 			if (DEBUG) for (RasterKey key : e.getValue()) {
-				LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.MEAT_OR_PASTURE), yieldRaster.get(key).getYieldMax(CropType.WHEAT)));
+				LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.PASTURE), yieldRaster.get(key).getYieldMax(CropType.WHEAT)));
 			}
 		}
 

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

@@ -11,7 +11,7 @@ public enum CommodityType {
 	PULSES("Pulses + (Total)", "pulses"),
 	SOYBEAN("Soyabeans", "soybean"),
 	STARCHY_ROOTS("Starchy Roots + (Total)", "starchyRoots"),
-	MEAT_OR_PASTURE("meatmilkeggs", "pasture_or_meat");
+	MEAT("meatmilkeggs", "meat");
 
 	private String faoName;
 	private String gamsName;

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

@@ -16,7 +16,7 @@ public enum CropType {
 	SOYBEAN("Soyabeans", "soybean"),
 	PULSES("Pulses + (Total)", "pulses"),
 	STARCHY_ROOTS("Starchy Roots + (Total)", "starchyRoots"),
-	MEAT_OR_PASTURE("meatmilkeggs", "pasture_or_meat");
+	PASTURE("meatmilkeggs", "pasture");
 
 	private String faoName;
 	private String gamsName;

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

@@ -30,7 +30,7 @@ public class YieldResponse {
 	
 	public double getFertParam() {
 		if (fertParm == 0)
-			fertParm = calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG));
+			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));
 		
 		return fertParm;
 	}