From ed37b6289e0cf8a555e5a173a9731d2c779040ec Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Wed, 8 Jul 2015 12:11:08 +0100
Subject: [PATCH] Split pasture and meat and fixes
---
GAMS/IntExtOpt.gms | 58 +++++++++++--------
src/ac/ed/lurg/country/CountryAgent.java | 2 +-
.../country/gams/GamsLocationOptimiser.java | 4 +-
.../lurg/country/gams/GamsLocationTest.java | 12 ++--
.../country/gams/GamsRasterOptimiser.java | 10 ++--
src/ac/ed/lurg/types/CommodityType.java | 2 +-
src/ac/ed/lurg/types/CropType.java | 2 +-
src/ac/ed/lurg/yield/YieldResponse.java | 2 +-
8 files changed, 50 insertions(+), 42 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 3aa6f79b..4ed20694 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -1,11 +1,11 @@
- 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;
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index d90227cd..4045f40d 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/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));
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 39a3bf30..0fc28992 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/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);
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationTest.java b/src/ac/ed/lurg/country/gams/GamsLocationTest.java
index 18d5dc88..d9112d76 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationTest.java
+++ b/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;
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index b47ff427..7d839add 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/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)));
}
}
diff --git a/src/ac/ed/lurg/types/CommodityType.java b/src/ac/ed/lurg/types/CommodityType.java
index aa51104e..ae9090a9 100644
--- a/src/ac/ed/lurg/types/CommodityType.java
+++ b/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;
diff --git a/src/ac/ed/lurg/types/CropType.java b/src/ac/ed/lurg/types/CropType.java
index 33543231..6f7d5a69 100644
--- a/src/ac/ed/lurg/types/CropType.java
+++ b/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;
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index 5aeb8d98..de53cad2 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/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;
}
--
GitLab