From c28b2aa72ba74d9cdd4ec71be5f20dce51894425 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Thu, 20 Aug 2015 11:15:32 +0100
Subject: [PATCH] Move seed and waste adjustment to demand
---
GAMS/IntExtOpt.gms | 16 +++--
config.properties | 14 ++++-
src/ac/ed/lurg/ModelConfig.java | 5 +-
src/ac/ed/lurg/ModelMain.java | 58 +++++++++----------
.../country/gams/GamsRasterOptimiser.java | 5 +-
src/ac/ed/lurg/demand/DemandManager.java | 2 +-
.../lurg/yield/LPJYieldResponseMapReader.java | 5 +-
7 files changed, 64 insertions(+), 41 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 500750e9..b02699c6 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -83,11 +83,7 @@ $gdxin
energy total input energy - energy;
POSITIVE VARIABLE area, fertI, irrigI, feedAmount, agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease;
-
- fertI.L(crop, location) = 0.5;
- irrigI.L(crop, location) = 0.5;
- area.L(crop, location) = previousArea(crop, location)
-
+
EQUATIONS
UNIT_ENERGY_EQ(crop, location) energy per area
YIELD_EQ(crop, location) yield given chosen intensity
@@ -174,6 +170,11 @@ $gdxin
* if fitToPreviousAreas is zero then we just optimise
if (fitToPreviousAreas eq 0,
+ fertI.L(crop, location) = 0.5;
+ irrigI.L(crop, location) = 0.5;
+ feedAmount.L(crop) = 0.0;
+ netImportAmount.L(import_crop) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
+ area.L(crop, location) = previousArea(crop, location)
SOLVE LAND_USE USING NLP MINIMIZING energy;
@@ -184,6 +185,11 @@ $gdxin
baseCropland = sum((location, crop_less_pasture), previousArea(crop_less_pasture, location));
while((count le 10 and not (abs(basePasture-newPasture) le 0.1 and abs(baseCropland-newCropland) le 0.1)),
+ fertI.L(crop, location) = 0.5;
+ irrigI.L(crop, location) = 0.5;
+ area.L(crop, location) = previousArea(crop, location);
+ feedAmount.L(crop) = 0.0;
+ netImportAmount.L(import_crop) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
SOLVE LAND_USE USING NLP MINIMIZING energy;
display pastureIncrease.l, pastureDecrease.l, cropDecrease.l, cropIncrease.l, unitEnergy.l;
diff --git a/config.properties b/config.properties
index a3348cf9..1bfc16aa 100644
--- a/config.properties
+++ b/config.properties
@@ -4,6 +4,16 @@ CLEANUP_GAMS_DIR=false
# Properties for testing
CHANGE_YIELD_DATA_YEAR=false
-#KEEP_DEMAND_FIXED=true
+KEEP_DEMAND_FIXED=true
+DEBUG_LIMIT_COUNTRIES=true
MAX_IMPORT_CHANGE=0.0
-SEED_AND_WASTE_FRACTION=0.1
\ No newline at end of file
+SEED_AND_WASTE_FRACTION=0.1
+
+
+DEBUG_LIMIT_COUNTRIES=true
+END_TIMESTEP=0
+TIMESTEP_SIZE=10
+POPULATION_AGGREG_LIMIT=60
+
+NUM_CEREAL_CATEGORIES=5
+NUM_PASTURE_CATEGORIES=1
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 1747fcac..3283c8da 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -139,6 +139,9 @@ public class ModelConfig {
public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 2.0);
public static final double TRANSPORT_LOSSES = getDoubleProperty("TRANSPORT_LOSSES", 0.2); // in international trade
-
+
+ public static final int NUM_CEREAL_CATEGORIES = getIntProperty("NUM_CEREAL_CATEGORIES", 5);
+ public static final int NUM_PASTURE_CATEGORIES = getIntProperty("NUM_PASTURE_CATEGORIES", 1);
+
public static final boolean DEBUG_LIMIT_COUNTRIES = getBooleanProperty("DEBUG_LIMIT_COUNTRIES", false);
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 2151c1e1..25764f8f 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -52,8 +52,8 @@ public class ModelMain {
private Map<CropType, Double> prevWorldPrices;
private RasterSet<IntensitiesItem> prevIntensityRaster;
private RasterSet<AreasItem> prevCropAreaRaster;
-
-
+
+
public static void main(String[] args) {
ModelMain theModel = new ModelMain();
theModel.setup();
@@ -63,11 +63,11 @@ public class ModelMain {
/* setup models, reading inputs, etc. */
private void setup() {
desiredProjection = RasterHeaderDetails.getGlobalHeaderFromCellSize(ModelConfig.CELL_SIZE_X, ModelConfig.CELL_SIZE_Y, "999");
-
+
BaseConsumpManager baseConsumpManager = new BaseConsumpManager();
compositeCountryManager = new CompositeCountryManager(baseConsumpManager);
demandManager = new DemandManager(ModelFitType.LOGISTIC, ModelConfig.SSP_SCENARIO, baseConsumpManager, compositeCountryManager);
-
+
countryBoundaryRaster = getCountryBoundaryRaster();
countryAgents = createCountryAgents(compositeCountryManager.getAll());
@@ -75,7 +75,7 @@ public class ModelMain {
prevWorldPrices = new HashMap<CropType, Double>();
for (CropType c : CropType.getImportedTypes())
prevWorldPrices.put(c, 0.3);
-
+
prevWorldPrices.put(CropType.STARCHY_ROOTS, 1.4);
}
@@ -95,7 +95,7 @@ public class ModelMain {
private void doTimestep(Timestep timestep) {
LogWriter.println(timestep.toString());
-
+
YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
// YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-50.0, -4.0);
@@ -107,15 +107,15 @@ public class ModelMain {
RasterSet<IntensitiesItem> globalIntensityRaster = new RasterSet<IntensitiesItem>(desiredProjection);
RasterSet<AreasItem> globalCropAreaRaster = new RasterSet<AreasItem>(desiredProjection);
RasterSet<IntegerRasterItem> globalLocationIdRaster = new RasterSet<IntegerRasterItem>(desiredProjection);
-
+
for (CountryAgent ca : countryAgents) {
LogWriter.println("Country " + ca.getCountry());
Collection<RasterKey> countryKeys = countryBoundaryRaster.getKeysFor(ca.getCountry());
YieldRaster countryYieldSurfaces = yieldSurfaces.getSubsetRasterForKeys(countryKeys);
-
+
GamsRasterOutput result = null;
-
+
// do the optimization
try {
result = ca.determineProduction(timestep, countryYieldSurfaces, prevWorldPrices);
@@ -135,7 +135,7 @@ public class ModelMain {
globalIntensityRaster.putAll(result.getIntensityRaster());
globalCropAreaRaster.putAll(result.getCropAreaRaster());
globalLocationIdRaster.putAll(result.getLocationIdRaster());
-
+
// Get values for world input costs
Map<CropType, CropUsageData> cropUsage = result.getCropUsageData();
@@ -149,7 +149,7 @@ public class ModelMain {
totalExportCommodities.incrementValue(c, -countryNetImports);
}
}
-
+
// Look at trade balance and adjust appropriately
for (CropType crop : CropType.getImportedTypes()) {
double prevPrice = prevWorldPrices.get(crop);
@@ -158,27 +158,27 @@ public class ModelMain {
LogWriter.printlnError(String.format("Price for %s not updated (still %.3f), as no imports or no exports for it", crop.getGamsName(), prevPrice));
continue;
}
-
+
double imports = totalImportCommodities.get(crop);
double exports = totalExportCommodities.get(crop) * (1.0-ModelConfig.TRANSPORT_LOSSES);
double adjustedPrice = updateMarketPrices(prevPrice, imports, exports);
LogWriter.println(String.format("Price for %s updated from %.3f (imports %.0f, exports %.0f) to %.3f ", crop.getGamsName(), prevPrice, imports, exports, adjustedPrice));
prevWorldPrices.put(crop, adjustedPrice);
}
-
-
+
+
// output results
outputTimestepResults(timestep, globalIntensityRaster, globalCropAreaRaster, globalLocationIdRaster, yieldSurfaces);
-
+
// keep last to allow interpolation
prevIntensityRaster = globalIntensityRaster;
prevCropAreaRaster = globalCropAreaRaster;
}
-
+
public double updateMarketPrices(double previousPrice, double demand, double supply) {
if (demand > 0 || supply > 0) {
double ratio;
-
+
if (demand > supply)
ratio = (demand-supply)/demand;
else
@@ -205,7 +205,7 @@ public class ModelMain {
return new BufferedWriter(fstream);
}
}
-
+
private double getTotalArea(LandCoverType lcType, RasterSet<AreasItem> globalCropAreaRaster) {
double totalArea = 0;
for (AreasItem item : globalCropAreaRaster.values()) {
@@ -219,13 +219,13 @@ public class ModelMain {
StringBuffer sbHeadings = new StringBuffer("Year, Cropland (Mha), Pasture (Mha), Natural (Mha)");
for (CropType crop : CropType.getImportedTypes())
sbHeadings.append(",Px_" + crop.getGamsName());
-
+
BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.TOTAL_LAND_COVER_FILE, sbHeadings.toString());
StringBuffer sbData = new StringBuffer();
sbData.append(String.format("%d,%.1f,%.1f,%.1f",
timestep.getYear(), getTotalArea(LandCoverType.CROPLAND, cropAreaRaster), getTotalArea(LandCoverType.PASTURE, cropAreaRaster), getTotalArea(LandCoverType.OTHER_NATURAL, cropAreaRaster)));
-
+
for (CropType crop : CropType.getImportedTypes() )
sbData.append(String.format(",%.3f", prevWorldPrices.get(crop)));
@@ -240,14 +240,14 @@ public class ModelMain {
private void outputTimestepResults(Timestep timestep, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster,
RasterSet<IntegerRasterItem> locationIdRaster, YieldRaster yieldSurfaces) {
-
+
writeMarketFile(timestep, cropAreaRaster);
-
+
if (ModelConfig.OUTPUT_FOR_LPJG) {
for (int outputYear : timestep.getYearsFromLast()) {
LogWriter.printlnError("Outputing Year: " + outputYear);
RasterSet<AreasItem> areasToOutput;
-
+
if (outputYear == timestep.getYear()) {
areasToOutput = cropAreaRaster;
}
@@ -261,12 +261,12 @@ public class ModelMain {
intermediateAreas.setup(prevCropAreaRaster, cropAreaRaster, timestep.getPreviousTimestep().getYear(), timestep.getYear(), outputYear);
areasToOutput = intermediateAreas;
}
-
+
LpjgOutputer lpjOutputer = new LpjgOutputer(outputYear, intensityRaster, areasToOutput, yieldSurfaces);
lpjOutputer.writeOutput();
}
}
-
+
outputLandCover(timestep.getYear(), cropAreaRaster, LandCoverType.CROPLAND);
outputLandCover(timestep.getYear(), cropAreaRaster, LandCoverType.PASTURE);
}
@@ -288,7 +288,7 @@ public class ModelMain {
CountryBoundaryRaster countryBoundaries = new CountryBoundaryRaster(desiredProjection, compositeCountryManager);
CountryBoundaryReader countryReader = new CountryBoundaryReader(countryBoundaries);
countryReader.getRasterDataFromFile(ModelConfig.COUNTRY_BOUNDARY_FILE);
-
+
return countryBoundaries;
}
@@ -298,16 +298,16 @@ public class ModelMain {
RasterSet<LandCoverItem> initLC = getInitialLandCover();
RasterSet<IrrigationCostItem> allIrrigationCosts = getIrrigationCosts();
Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap = new CropUsageReader(compositeCountryManager).getCommodityData();
-
+
for (CompositeCountry cc : countryGrouping) {
-
+
// DEBUG code
if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {
if (!(cc.getName().equals("United States of America") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
continue;
}
}
-
+
List<RasterKey> keys = countryBoundaryRaster.getKeysFor(cc);
RasterSet<LandCoverItem> initCountryLC = initLC.popSubsetForKeys(new RasterSet<LandCoverItem>(initLC.getHeaderDetails()), keys);
RasterSet<IrrigationCostItem> irrigationCosts = allIrrigationCosts.popSubsetForKeys(new RasterSet<IrrigationCostItem>(allIrrigationCosts.getHeaderDetails()), keys);
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 906d62e4..aa4369e4 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -8,6 +8,7 @@ import java.util.Map;
import java.util.Map.Entry;
import java.util.Set;
+import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.landuse.AreasItem;
import ac.ed.lurg.landuse.IntensitiesItem;
import ac.ed.lurg.landuse.IrrigationCostItem;
@@ -255,8 +256,8 @@ public class GamsRasterOptimiser {
}
- int numCerealCats = 5;
- int numPastureCats = 1;
+ int numCerealCats = ModelConfig.NUM_CEREAL_CATEGORIES;
+ int numPastureCats = ModelConfig.NUM_PASTURE_CATEGORIES;
int thisShouldLookAtCropsOtherThanJustWheat; // need to consider other crops, and perhaps other yieldTypes as well
List<Double> wheatlDivisions = getDivisions(yieldRaster, CropType.WHEAT, numCerealCats);
diff --git a/src/ac/ed/lurg/demand/DemandManager.java b/src/ac/ed/lurg/demand/DemandManager.java
index 772c0c11..e666003d 100644
--- a/src/ac/ed/lurg/demand/DemandManager.java
+++ b/src/ac/ed/lurg/demand/DemandManager.java
@@ -58,7 +58,7 @@ public class DemandManager {
double bioenergy = getBioenergyDemand(c, year, dc.getCommodityType());
// LogWriter.println(String.format("Country %s comm %s: %f", country, dc.getCommodityType(), bioenergy));
- double newDemand = d + bioenergy;
+ double newDemand = (d + bioenergy) * (1.0 + ModelConfig.SEED_AND_WASTE_FRACTION);
Double prevDemand = demandMap.get(dc.getCommodityType());
if (prevDemand != null)
diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
index dbe550b6..74c14780 100644
--- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
+++ b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
@@ -72,7 +72,7 @@ public class LPJYieldResponseMapReader extends AbstractTabularRasterReader<Yield
}
else {
// TeWWhi TeWWme TeWWlo TeSWhi TeSWme TeSWlo TeCohi TeCome TeColo TrRihi TrRime TrRilo
- double adjFactor = 10 * (1.0 - ModelConfig.SEED_AND_WASTE_FRACTION) ; // 10 for kg/m2 to t/ha
+ double adjFactor = 10; // 10 for kg/m2 to t/ha
boolean isSpringWheat = (getValueForCol(rowValues, "teSWhi") > getValueForCol(rowValues, "TeWWhi"));
item.setIsSpringWheat(isSpringWheat);
@@ -81,6 +81,9 @@ public class LPJYieldResponseMapReader extends AbstractTabularRasterReader<Yield
for (IrrigationRate irrig : IrrigationRate.values()) {
YieldType yieldType = YieldType.getYieldType(fert, irrig);
+ if (fert.equals(FertiliserRate.NO_FERT))
+ adjFactor *= 0.5;
+
String fertIrrigString = irrig.getId() + fert.getId();
double ww = getValueForCol(rowValues, "TeWW" + fertIrrigString) * adjFactor;
double sw = getValueForCol(rowValues, "TeSW" + fertIrrigString) * adjFactor;
--
GitLab