diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 603cdfc02b9cec9a426365f15694527dc4444f42..f73e4f799efe6442804473cac31d1a4922f2cb39 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -40,12 +40,17 @@ $gdxin
SCALAR minDemandPerCereal / 0.1 /
- 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 / ;
+ 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
+ soybean 0.7
+ pasture 1 / ;
+* pasture yield is done in DM terms
+
+
+ PARAMETER zeroIntensityYieldFrac(crop) fraction of no fertiliser and irrigiation yield with zero intensity. Most or all arable crops this is 0 but not for pasture
+ / pasture 0.6 / ;
VARIABLES
area(crop, location) total area for each crop - Mha
@@ -53,7 +58,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,7 +80,7 @@ $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
+ 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
@@ -89,8 +94,8 @@ $gdxin
ENERGY_EQ total energy objective function;
UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= fertI(crop, location) +
- ((1+irrigI(crop, location) ** 2) * irrigCost(location)) +
- (1+otherIntensity(crop, location) ** 2);
+ ((1+irrigI(crop, location) ** 2) * irrigCost(location)) -1 +
+ (1+otherIntensity(crop, location) ** 2) - 1;
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yieldNone(crop, location) +
@@ -98,14 +103,14 @@ $gdxin
(yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) +
(yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
((fertI(crop, location) + delta) ** fertParam(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location))
- ) * otherIntensity(crop, location);
+ ) * (zeroIntensityYieldFrac(crop) + otherIntensity(crop, location) * (1-zeroIntensityYieldFrac(crop)));
* YIELD_EQ(crop, location) .. yield(crop, location) =E= (
* yieldNone(crop, location) +
* (yieldFertOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-10*fertI(crop, location)) +
* (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-5*irrigI(crop, location)) +
* (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);
+* ) * (zeroIntensityYieldFrac(crop) + otherIntensity(crop, location) * (1-zeroIntensityYieldFrac(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);
@@ -115,12 +120,12 @@ $gdxin
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') * sum(location, area('pasture', location) * yield('pasture', location)) =G=
+ 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;
-* MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) .. otherIntensity(crop, location) =L= 10;
+ MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) .. otherIntensity(crop, location) =L= 1;
TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop, area(crop, location));
@@ -147,8 +152,8 @@ $gdxin
MODEL LAND_USE /ALL/ ;
SOLVE LAND_USE USING NLP MINIMIZING energy;
-display feedAmount.l, minFeedRate, demand, netImportAmount.l;
-
+display meatEfficency, feedDM, feedAmount.l, area.l, yield.l, demand, netImportAmount.l;
+
Scalar ms 'model status', ss 'solve status';
ms=LAND_USE.modelstat;
ss=LAND_USE.solvestat;
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 49ab8c57ce5b41dd7383bcba968289784b653748..50b8487bbe2ed898ba4a5ab4d5ed0a8270342a8f 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -163,7 +163,7 @@ public class ModelMain {
RasterSet<IrrigationCostItem> allIrrigationCosts = getIrrigationCosts();
Map<Country, Map<CropType, CropUsageData>> cropUsageDataMap = CropUsageData.readCommodityData();
- HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("French Polynesia", "Cabo Verde", "Samoa", "Saint Vincent and the Grenadines"));
+ // HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("French Polynesia", "Cabo Verde", "Samoa", "Saint Vincent and the Grenadines"));
for (Map.Entry<Country, List<RasterKey>> entry : countryToKeysMap.entrySet()) {
Country country = entry.getKey();
@@ -174,10 +174,8 @@ public class ModelMain {
if (!country.getCountryName().equals("United States of America")) {
continue;
}
-
-
- if (countryExclusionList.contains(country.getCountryName())) {
+ if (demandManager.getPopulation(country, 2010) < 5) { // countryExclusionList.contains(country.getCountryName()
LogWriter.printlnError("Skipping " + country);
continue;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index c5a2f321a81759a7d8b668bb7d18d29c56c4fea2..abe96dd422bc76e3a8fa96afda0b60a1aa677f94 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -58,15 +58,15 @@ public class GamsCountryInput {
// }
public double getMeatEfficiency() {
- return 1; // this is already handled by the feed conversion efficiency for each animal product
+ return 0.5; // this is already handled by the feed conversion efficiency for each animal product
}
public double getMaxLandUseChange() {
- return 0.05;
+ return 0.95;
}
public double getLandChangeEnergy() {
- return 0.1;
+ return 0.001;
}
public double getMinFeedRate() {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 6b82ca3bdb41f571d9897d167202ee2cae0c0b41..1daa46b61faee818d948c0cfa8a2de7069272829 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -85,6 +85,8 @@ public class GamsLocationOptimiser {
GAMSParameter prevCropP = inDB.addParameter("previousArea", 2);
GAMSParameter landP = inDB.addParameter("suitableLandArea", 1);
+ double totalAgriLand = 0;
+
for (Map.Entry<Integer, ? extends AreasItem> entry : inputData.getPreviousAreas().entrySet()) {
Integer locationId = entry.getKey();
String locString = Integer.toString(locationId);
@@ -97,13 +99,15 @@ public class GamsLocationOptimiser {
v.add(cropType.getGamsName());
v.add(locString);
prevCropP.addRecord(v).setValue(d);
+ totalAgriLand += d;
}
double suitableLand = areasItem.getSuitableLand();
if (DEBUG) LogWriter.println(String.format(" %d %15s,\t %.1f", locationId, "suitableLand", suitableLand));
landP.addRecord(locString).setValue(suitableLand);
}
-
+ if (DEBUG) LogWriter.println(String.format(" Total agricultural\t %.1f", totalAgriLand));
+
if (DEBUG) LogWriter.println("\nIrrigation cost");
GAMSParameter irrigCostP = inDB.addParameter("irrigCost", 1);
for (Entry<Integer, ? extends IrrigationCostItem> entry : inputData.getIrrigationCosts().entrySet()) {
@@ -118,7 +122,7 @@ public class GamsLocationOptimiser {
GamsCountryInput countryInput = inputData.getCountryInput();
addCommodityMapParm(inDB.addParameter("demand", 1), countryInput.getProjectedDemand());
- if (DEBUG) LogWriter.println("\nYield (fert/irrig) None/None, Max/None, None/Max, Max/Max, [fert p], [irrig p]");
+ if (DEBUG) LogWriter.println("\nYield (fert/irrig) None/None, Max/None, None/Max, Max/Max,\t [fert p],\t [irrig p]");
GAMSParameter yNoneP = inDB.addParameter("yieldNone", 2);
GAMSParameter y_fert = inDB.addParameter("yieldFertOnly", 2);
GAMSParameter y_irrig = inDB.addParameter("yieldIrrigOnly", 2);
@@ -132,8 +136,8 @@ public class GamsLocationOptimiser {
YieldResponsesItem yresp = entry.getValue();
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)));
+ if (DEBUG) LogWriter.println(String.format("%d %15s,\t %.1f,\t %.1f, \t %.1f,\t %.1f,\t\t [%.2f],\t [%.2f]",
+ locationId, crop.getGamsName(), yresp.getYieldNone(crop), yresp.getYieldFertOnly(crop), yresp.getYieldIrrigOnly(crop), yresp.getYieldMax(crop), yresp.getFertParam(crop), yresp.getIrrigParam(crop)));
Vector<String> v = new Vector<String>();
v.add(crop.getGamsName());
@@ -235,10 +239,16 @@ public class GamsLocationOptimiser {
}
prevCropType = cropType;
}
- if (DEBUG) LogWriter.println(String.format("\nTotal area= %.1f", totalArea));
-
+
+ netImport = varNetImports.findRecord(CropType.MEAT.getGamsName()).getLevel();
+ cropUsageData.put(CropType.MEAT, new CropUsageData(0.0, netImport));
+ if (DEBUG) {
+ LogWriter.println(String.format("\n%s:\t\t\t\t\tnetImports= %.3f", CropType.MEAT.getGamsName(), netImport));
+ LogWriter.println(String.format("\nTotal area= %.1f", totalArea));
+ }
+
GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, cropUsageData);
- return results ;
+ return results;
}
private void addScalar(GAMSParameter param, double val) {
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index c05d29b6350ce0fae9069c846032a57758fc898b..31096c523521b6ee115871cb6ce3b24ca309aedd 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -67,7 +67,7 @@ public class GamsRasterOptimiser {
// if (DEBUG) LogWriter.println("Processing location id " + locId);
double shortfall = allocAllAreasForAgg(newAreaRaster, prevAreaRaster, keys, cropChangeTotals);
- if (shortfall > 0)
+ if (shortfall > 0.00001)
LogWriter.printlnError("This should never happen, due to GAMS constraint. Not able to incorporate all changes, as not enough forest or natural areas left: " + shortfall);
}
@@ -119,7 +119,7 @@ public class GamsRasterOptimiser {
ratioOfChanges = (prevForest + prevNatural) / netAddedCrop;
forestRemoved = prevForest;
naturalRemoved = prevNatural;
- LogWriter.println("Not able to incorporate all changes for raster location: " + netAddedCrop);
+ if (DEBUG) LogWriter.println("Not able to incorporate all changes for raster location: " + netAddedCrop);
}
else if (prevForest < forestRemoved) {
forestRemoved = prevForest;
diff --git a/src/ac/ed/lurg/demand/DemandManager.java b/src/ac/ed/lurg/demand/DemandManager.java
index ace23a3aa4551b84c42c10b4bec6374b1d24446e..061e2ae48141602b274549a930031b88b8cfec05 100644
--- a/src/ac/ed/lurg/demand/DemandManager.java
+++ b/src/ac/ed/lurg/demand/DemandManager.java
@@ -47,6 +47,16 @@ public class DemandManager {
}
return demandMap;
}
+
+ public double getPopulation (Country country, int year) {
+ SspData sd = sspManager.get(ssp_scenario, year, country);
+ if (sd == null) {
+ LogWriter.printlnError("No ssp data for " + ssp_scenario + ", " + year + ", " + country.getCountryName());
+ return 0;
+ }
+ else
+ return sd.getPopulation();
+ }
public String getModelAndScenarioDescription() {
return fitType + " " + ssp_scenario;
diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
index 9197fc6883b402f55d05a8c35ca5df4fb5ef500d..cee07b613ad023410a95f14b25650a24a5845dc2 100644
--- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
+++ b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
@@ -123,8 +123,8 @@ public class LPJYieldResponseMapReader {
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);
+ item.setYield(noIrrigYieldType, CropType.PASTURE, record.teSW * 2);
+ item.setYield(maxIrrigYieldType, CropType.PASTURE, record.teSW * 2.1);
}
class YieldRecord {