From bb4a660e882e90fae5c9a6ffea189a40e03363b2 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Wed, 8 Jul 2015 20:56:39 +0100
Subject: [PATCH] Fixes, including handling pasture and meat separately
---
GAMS/IntExtOpt.gms | 46 ++++++++++---------
src/ac/ed/lurg/ModelConfig.java | 2 +-
.../country/gams/GamsLocationOptimiser.java | 9 ++--
.../country/gams/GamsRasterOptimiser.java | 6 +--
src/ac/ed/lurg/types/CropType.java | 40 ++++++++++++++--
.../lurg/yield/LPJYieldResponseMapReader.java | 41 +++++++++++------
src/ac/ed/lurg/yield/YieldResponse.java | 7 ++-
src/ac/ed/lurg/yield/YieldResponsesItem.java | 4 ++
8 files changed, 106 insertions(+), 49 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 4ed20694..603cdfc0 100644
--- a/GAMS/IntExtOpt.gms
+++ b/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;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 1d999e92..63a07041 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/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";
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 0fc28992..6b82ca3b 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/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));
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 7d839add..c05d29b6 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/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));
}
diff --git a/src/ac/ed/lurg/types/CropType.java b/src/ac/ed/lurg/types/CropType.java
index 6f7d5a69..8ed5b07e 100644
--- a/src/ac/ed/lurg/types/CropType.java
+++ b/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;
+ }
}
diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
index 343b7abe..9197fc68 100644
--- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
+++ b/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 {
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index de53cad2..ad89e96c 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/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;
}
diff --git a/src/ac/ed/lurg/yield/YieldResponsesItem.java b/src/ac/ed/lurg/yield/YieldResponsesItem.java
index 6a068ca7..cc8f6195 100644
--- a/src/ac/ed/lurg/yield/YieldResponsesItem.java
+++ b/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;
}
--
GitLab